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SECTION  I 


INTRODUCTION 


Predictions  of  gas  turbine  combustor  performance  and  pollutant  emission 
characteristics  require  modeling  procedures  possessing  a high  degree  of 
sophistication.  Past  attempts  at  modeling  combustion  systems  have  been 
largely  frustrated  by  the  complexity  of  the  coupled  hydrodynamic  and  chemical 
processes.  The  difficulty  can  be  largely  attributed  to  the  lack  of  under- 
standing of  the  flow  processes  which,  through  the  exchange  of  heat,  mass  and 
momentum,  can  directly  relate  to  pollutant  formation  and  combustion  efficien- 
cy. For  example,  swirling  flow  has  been  shown  to  have  an  important  influence 
on  the  stability  and  combustion  intensity  of  flames  (Ref.  1)  as  well  as  on 
residence-time  distributions  (Ref.  2)  in  combustors  which,  in  turn,  can  be 
related  to  combustor  performance  and  efficiency  as  well  as  to  emission  charac- 
teristics (Refs.  3 and  4).  Techniques  employed  in  modeling  combustor  flow 
processes  have  generally  been  highly  si  olified,  particularly  flow  modeling 
techniques  where  stirred  reactor  concept  and  one -dimensional  assumptions  are 


1.  Beer,  J.  M.  and  N.  A.  Chigier:  Stability  and  Combustion  Intensity  of 

Pulverized  Coal  Flames  - Effect  of  Swirl  and  Impingement.  Journal  of  the 
Institute  of  Fuel,  December  1969- 

2.  Beer,  J.  M.  and  W.  Leucker:  Turbulent  Flames  in  Rotating  Flow  Systems. 

Paper  No.  Inst.  F-NAFTC-7,  North  American  Fuel  Technology  Conference, 
Ottawa,  Canada . 1970. 

3.  Beer,  J.  M.  and  J.  B.  Lee:  The  Effects  of  Residence  Time  Distribution  on 

the  Performance  and  Efficiency  of  Combustors.  The  Combustion  Institute, 
1965,  pp.  1187-1202. 

k.  Marteney,  P.  J.:  Analytical  Study  of  the  Kinetics  of  Formation  of  Nitrogen 

Oxide  in  Hydrocarbon  - Air  Combustion.  Combustion  Science  and  Technology. 
Vol.  1,  1970,  pp.  37-^5 . 
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employed  (Refs.  5 through  10).  Chemistry  is  frequently  modeled  by  assuming 
equilibrium  hydrocarbon  fuel  decomposition  and  two  phase  flow  effec* s are 
seldom  considered.  In  some  more  recent  modeling  attempts,  for  example. 
Fletcher  and  Heywood  (Ref.  5)  and  Hammond  and  Mellor  (Refs.  6 and  7),  the 
stirred  reactor  concept  is  employed  to  assess  the  effect  of  residence-time  on 
combustion  behavior  and  to  predict  pollutant  emissions  in  gas  turbines. 

Droplet  vaporization  and  burning  were  neglected  in  these  studies;  however,  a 
quasi -global  finite-rate  hydrocarbon  combustion  mechanism  was  employed  by 
Hammond  and  Mellor  to  model  the  chemistry.  In  related  work,  Roberts,  et  al., 
(Ref.  8),  in  an  attempt  to  predict  nitrogen  oxide  formation  in  gas  turbine 
combustors,  subdivided  the  combustor  into  three  regions;  one  corresponding 
to  the  central  recirculation  portion  of  the  upstream  zone;  a second  repre- 
senting the  flow  region  surrounding  the  recirculation  zone  which  was  inter- 
preted to  be  a one -dimensional  reacting  zone;  and  the  third  downstream  zone 
modeled  as  a one-dimensional  region.  Both  finite-rate  and  equilibrium  hydro- 
carbon chemistry  models  were  considered.  It  is  interesting  that  little  dif- 
ference in  the  predicted  NO  levels  was  noted  in  their  results  between  the 
equilibrium  and  finite-rate  hydrocarbon  cases.  A more  recent  analysis  di- 
rected toward  low  power  application  by  Mosier,  et  al.,  (Ref.  9)  basically 
extended  the  work  of  Ref.  8 through  the  use  of  a more  sophisticated  finite- 
rate  hydrocarbon  chemistry  model  obtaining  trends  in  agreement  with  experi- 
mental data.  The  modular  approach  proposed  by  Edelman  and  Economos  (Ref.  10) 
is  an  attempt  at  formulating  a general  analytical  procedure  for  predicting 
combustor  behavior  by  treating  the  various  critical  combustor  processes  on  an 
individual  basis  or  coupled  as  a function  of  operating  conditions.  Difficulty 
with  the  approach  lies  with  its  method  of  accounting  for  recirculation  (a 
stirred  reactor  is  presently  used)  and  its  inability  to  provide  a unified 
description  of  the  burner  under  a given  set  of  operating  conditions. 


5.  Fletcher,  R.  S.  and  J.  B.  Heywood:  A Model  for  Nitric  Oxide  Emission  from 

Aircraft  Gas  Turbine  Engines.  AIAA  Paper  71-123,  1971. 

6.  Hammond,  D.  C..  Jr.  and  A.  M.  Mellor;  Analytical  Predictions  of  Emissions 
from  and  Within  an  Allison  J-33  Combustor.  Combustion  Science  and  Tech- 
nology, Vol.  6,  1973,  pp.  279-286. 

7.  Hammond,  D.  C.,  Jr.  and  A.  M.  Mellor:  Analytical  Calculations  fbr  the 

Performance  and  Pollutant  Emissions  of  Gas  Turbine  Combustors.  Combustion 
Science  and  Technology.  Vol.  4,  1971,  pp.  101-112. 

1‘  • Roberts,  R. , L.  D.  Aceto . R.  Keilback,  D.  P.  Teixeira,  and  J.  M.  Bonnell: 

An  Analytical  Model  for  Nitric  Oxide  Formation  in  a Gas  Turbine  Combustion 
Chamber.  AIAA  Paper  No.  7I-715,  1971. 

Q.  Mosier,  S.  A.,  R.  Roberts,  and  R.  E.  Henderson:  Development  and  Verifi- 

cation of  an  Analytical  Model  for  Predicting  Emissions  from  Gas  Turbine 
Engine  Combustors  During  Low  Power  Operation.  4lst  Meeting  Propulsion  and 
Energetics  Panel  of  AGARD,  1973* 

10.  Edelman,  R.  and  C.  Economos:  A Mathematical  Model  for  Jet  Engine  Combustor 

Pollutant  Emissions.  AIAA  FAper  No.  71-714,  1971. 
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The  foregoing  methods  are  lacking  primarily  in  their  ability  to  proper 
account  for  mixing  phenomena  occurring  in  the  reverse  flow  recirculation  con- 
of  combustion  devices.  It  has  recently  become  feasible,  however,  to  treat 
more  rigorously  flows  having  recirculation  zones,  by  numerical  solution  <•! 
elliptic  equations  governing  combusting  flows.  For  example,  one  such  nui.v  ; h-a. 
method  based  on  an  explicit  point -by -point  relaxation  procedure  has  bt--n  sug- 
gested by  Gosman,  et  al.,  (Ref.  11).  Sample  calculations  of  a representative 
gas  turbine  combustor  flow  have  been  computed  at  United  Technologies  Research 
Center  using  the  Gosman,  et  al.  method,  to  demonstrate  the  feasibility  of 
making  computations  in  the  recirculating  zones  of  combustion  chamber:  . Th 
results  obtained  with  this  procedure  demonstrated  qualitative  agreement  with 
experimental  observations  (Ref.  12). 

Although  these  results  were  very  encouraging,  the  slow  convergence 
properties  of  the  Gosman,  et  al . procedure,  arising  primarily  from  the  use  o: 
a point -by-point,  relaxation  technique,  became  apparent.  Consequently,  an 
improved  numerical  procedure  was  developed  at  UTRC  for  solving  combusting 
flows  containing  recirculation  zones.  The  UTRC  procedure  is  an  implicit  com- 
putational scheme,  and  is  novel  in  that  residuals  are  relaxed  simultaneously 
throughout  the  entire  flow  field,  rather  than  one  at  a time,  as  is  character- 
istic of  the  explicit  point  methods.  Under  a joint  AFAPL/FAA  contract,  the 
UTRC  Field  R- luxation  Elliptic  Procedure  (FRET')  ba. ed  on  the  rigorous  s> lution 
of  the  governing  equations  was  further  developed  and  used  to  predict  the  per- 
formance and  emission  characteristics  of  can-annular  and  annular  gas  turbine 
combustors  (Ref.  13)-  Further  development  of  the  procedure  and  the  physical 
models  is  presently  being  carried  out  under  EPA  Contract  No.  68-02-1873.  Sir- 
UTRC  method  solves  the  axisymmetric  time -averaged  Navier-Stokes  equations 
including  the  effects  of  turbulence,  chemistry,  radiation,  and  droplet  vapor- 
ization. The  FREP  code  has  given  reasonable  predictions  for  combustor  flow 
fields  which  have  axial  symmetry;  however,  significant  circumferential  asym- 
metry is  present  in  many  aircraft  combustor  flow  fields  and  a realistic 
approach  must  consider  this  complication.  Such  a general  approach  involves 
solution  of  the  time  averaged  three-dimensional  Navier-Stokes  equations  in- 
cluding turbulence,  chemical  reactions,  radiation  transport  and  droplet  vapor- 
ization. Recently,  for  instance,  Patankar  and  Spalding  (Ref.  l41  have  develop- 1 


11.  Gosman,  A.  D.,  W.  M.  Pun,  A.  K.  Runchal,  P.  R.  Spalding,  and  M.  Wolfsht  <? if. : 
Heat  and  Mass  Transfer  in  Recirculating  Flows.  Academic  Press,  N°w  York, 
1969. 

12.  Anasoulis,  R.  F. : Computations  of  the  Flow  in  a Combustor.  United 

Aircraft  Research  laboratories  Report  Kll.0885-1 . November  ! >71. 

1.3*  Anasoulis,  R.  F.  , H.  McDonald  and  R.  C.  Buggeln;  Development  of  a 
combustor  Flow  Analysis.  part  I:  Theoretical  'tudi-'- . Air  Force 

t’ropulsion  laboratory  Report  AFAPI.-TR -71-98,  Part  , January  I >7... 

14 . Patankar,  S.  V.  and  D.  P.  Spalding:  A Computer  Model  for  hree-'  i tn- -i . . . 

Flow  in  Furnaces.  Fourteenth  Symposium  (International  on  v>mt  • ; 'ion.  ihe 
Combustion  Ins?  itut  . 1971.  pp.  > . o-t>l4. 
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a line  relaxation  procedure  for  computation  of  steady  three-dimensional 
combusting  flows  in  cartesian  geometries;  however,  this  procedure  is  not 
generally  available  and  few  details  of  this  technique  are  known. 


■ 


The  objective  of  the  present  investigation  was  to  extend  an  existing  and 
relatively  efficient  UTRC  three-dimensional  Navier-Stokes  calculation  proce- 
dure (Ref.  15)  so  that  it  would  be  able  to  compute  combusting  flows.  In 
particular,  the  procedure  developed  includes  a simple  mixing  length  turbu- 
lence model,  a pseudo-kinetic  hydrocarbon  chemistry  model,  a liquid  droplet 
vaporization  model,  a single  frequency  radiation  model  and  a finite  rate 
nitrogen  oxide  chemistry  model.  The  intent  in  the  development  of  the  flow 
models  described  above  was  to  eliminate  undue  complexity  and  sophistication, 
simultaneously  providing  a reasonably  good  framework  within  which  refinements 
could  be  easily  implemented  at  a future  time,  if  warranted  by  comparisons 
with  experimental  data. 


15 - Briley,  W.  R.  and  H.  McDonald:  An  Implicit  Numerical  Method  for  the 

Multidimensional  Compressible  Navier-Stokes  Equations.  United  Aircraft 
Research  laboratories  Report  M9H 363-6,  November  1973. 
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SECTION  II 


TSEORETICAL  ANALYSIS 


Approach 


The  flow  regime  considered  in  the  present  study  is  a steady  or  unsteady 
gas-phase  flow  with  hydrocarbon  chemistry,  droplet  vaporization  and  burnin,', 
and  radiation  transport.  An  equilibrium  mixing  length  model  with  a general- 
ized eddy  viscosity  is  used  to  specify  the  turbulent  momentum  i luxes  (Reynolds’ 
stress)  in  the  time-averaged  Navier-Stokes  equations.  The  turbulent  fluxes  of 
enthalpy  and  chemical  species  are  determined  by  specifying,  turbulent  exchange 
coefficients  using  values  for  effective  Prandtl  and  Schmidt  numbers  taken  from 
knowledge  of  turbulent  flow  of  gases  and  gas  mixtures.  In  addition,  a chemis- 
try model,  a droplet  model,  and  a radiation  model  are  necessary  to  include  the 
effects  of  chemical  reactions,  two  phase  flow,  and  radiation  emission  and 
absorption  on  the  combustor  flow  field. 


A computational  method  is  required  to  solve  the  complex  system  of 
equations  obtained  for  combusting  flows  using  the  physical  models  describ'd 
above.  The  computational  procedure  must  be  capable  of  treating  the  flow  fielu 
resulting  from  the  mixing  and  chemical  reaction  of  the  appropriate  chemical 
species,  and  from  large  gradients  in  flow  properties  caused  by  the  combustion 
process.  The  Multidimensional  Implicit  Nonlinear  Time-dependent  (MINT)  tech- 
nique developed  by  Briley  and  McDonald  (Ref.  15)  for  the  compressible  Navier- 
Stokes  equations  is  well  suited  for  application  to  the  complex  equations 
governing  combusting  flows  and  was  employed  in  the  present  study.  The  MINT 
procedure  is  particularly  attractive  in  view  of  its  economic  computer  stora;  *- 
requirements  (only  a portion  of  the  flow  is  required  in  core  at  any  given 
time)  and  its  high  efficiency  relative  to  other  available  schemes.  The  re- 
sulting computer  code  is  used  to  compute  time-mean-average  velocities,  tem- 
peratures, pressures  and  species  concentrations  within  a selected  combustor 
with  discrete  inlet  injection  ports,  and  the  results  are  compared  with  the 
available  experimental  measurements  to  allow  an  evaluation  of  the  procedure 
and  the  analytical  modeling  techniques. 


Governing  Equations 

Under  consideration  is  the  flow  of  a turbulent  chemically  reactin'  mul ti - 
component  mixture  with  heat  and  mass  transport.  The  governing  system  of  par- 
tial differential  equations  describing  the  combustion  process  is  based  on  the 
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conservation  laws  of  mass,  momentum,  energy,  and  chemical  species  (Ref.  16). 

For  simplicity  these  equations  are  expressed  in  vector  notation  oelow  and  all 
quantities  are  nond imens ional . Velocities  are  normalized  by  Up,  density  by 
Dp,  enthalpy  by  hp,  temperature  by  Tp,  molecular  weights  by  Wp,  pressure  by 
Pp  = PpRgTpZp(Zp  - 107Wp),  dynamic  viscosity  by  (j,p,  radiation  energy  flux  by 
qpp,  and  time  by  (L/Up)  where  L is  the  reference  length.  Coupling  between 
concentration  and  thermal  gradients  (Soret  and  Dufour  effects),  pressure  gra- 
dient diffusion,  body  forces  and  bulk  viscosity  are  all  assumed  to  be  negligi- 
ble. In  addition,  Fick's  law  is  presumed  valid  which  implies  equal  binary 
diffusion  coefficients  for  each  pair  of  species  in  the  mixture  (see,  e.r.,  Ref. 
17).  The  resulting  set  of  time-averaged  equations  is 


Continuity 


= - V (P  u) 


Conservation  of  (las  Phase  Species 


Conservation  of  Liquid  Phase  Species 


* -V  (p3f,)+  j-  v.ffpVfj)  + Rj 


Conservation  of  Momentum 


^a)..7gpuC)-lL  Jp.itM-fii’'  Km 


(V  0)1  (<0 


Conservation  of  Energy 


^h!.-7.(p0h)+-£e_  |e_  ,-Lv  (rhVH) 

+ Re  7 [(rm_rh)  I h,Vm,]-  d-  V - g 
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The  mean  flow  rate  of  strain  tensor  in  Eq.  (4)  is  given  by 


e:  ±.  [ ( Vu ) +( Vu)T] 


(6) 


The  necessary  thermodynamic  relationships  are 


P=  PTZ 
Wm  f W, 


HrZ  m-  h:  + -L  ~r~  tu-u) 

i 2 hn 


(7a) 

(7b) 

(8) 


and  the  enthalpy  of  species  i is 


hi  r/  cp.(T)  dT 


f 


(9) 


Note  that  the  heat  of  formation  (in  ) of  species  i does  not  appear  in  the 
definition  Eq.  (9) 5 but  rather  has  been  included  explicitly  in  the  energy  con- 
servation equation  (5).  This  formulation  was  found  beneficial  in  the  present 
procedure  in  order  to  reduce  numerical  difficulties  due  to  truncation  errors. 
The  kinetic  heating  terms  in  the  energy  equation  which  are  not  significant  for 
the  Mach  numbers  under  consideration,  have  been  neglected. 


In  order  to  solve  the  above  system  of  equations,  in  addition  to  boundary 
conditions,  it  is  necessary  to  specify  the  turbulent  exchange  coefficients 
Heffj  r^,  r , T ; the  rate  of  production  of  species  i due  to  chemical,  reaction, 
r • , the  source  due  to  vaporization  of  liquid  droplets  from  particle  class  j, 


the  droplet  source  term  Rf;  and  the  radiation  energy  source  term, 

In  the  present  analysis  since  effective  Prandtl  and  Schmidt  numbers 


1 : 

s i»_j» 

“V’qR* 

are  defined  from  knowledge  of  turbulent  flows  of  gases  and  gas  mixtures,  only 
the  turbulent  momentum  exchange  coefficient,  must  be  specified.  The 

energy  and  species  exchange  coefficients  are  obtained  from  the  relations 


rh  = 


Meff 

preff 

(10) 

_ Meff 

(ID 

Sceff 

A turbulence  model  is  employed  to  define  the  effective  viscosity,  y, 

Similarly,  a chemistry  model  is  employed  to  specify  the  production  rate  r-j  for 
hydrocarbon  and  nitrogen  oxide  chemistry,  a droplet  vaporization  model  is 
employed  for  the  source  terms  R.  and  s-  . , and  a radiation  transport  model 

_ *.  J j J 

serves  to  specify  v • q^ . These  models  will  be  discussed  in  detail  subsequently. 
Geometry  and  Coordinate  System 

The  governing  vector  equations  presented  above  must  be  written  in  a 
coordinate  system  appropriate  for  combustor  flow.  In  the  present  study  con- 
sideration was  directed  primarily  toward  three-dimensional  flown  in  axially 
symmetric  combustor  geometries  with  a discrete  circumferential  distribution 
of  air  and  fuel  injection  ports.  Rectangular  duct  geometries  may  also  be 
treated  quite  easily  within  this  framework.  To  obtain  axisymmetric  coordi- 
nates of  sufficient  generality,  a two-dimensional  orthogonal  curvilinear  co- 
ordinate system  is  rotated  about  an  axis  of  symmetry  to  produce  the  desired 
geometry  (Fig.  l).  The  axisymmetric  coordinates  are  derived  from  a general 
system  of  orthogonal  curvilinear  coordinates  x-^,  x2,  x^  with  metric  coeffi- 
cient hj,  h2,  h^ . The  vector  operations  necessary  for  deriving  the  governing 
differential  equations  in  this  coordinate  system  from  the  vector  equations  (1) 
through  (5)  may  be  found  in  Ref.  18.  These  are  summarized  in  Appendix  A for 
completeness.  These  vector  relations  may  be  substituted  into  equations  (1) 
to  (5)  to  obtain  the  governing  equations  in  the  generalized  coordinate  sys- 
tem. The  system  of  axisymmetric  coordinates  are  obtained  by  taking  h2  = r 
and  x2  = 9,  where  r is  the  radial  distance  from  the  axis  of  symmetry  and  9 
is  the  angular  coordinate  (Fig.  l). 

Turbulence  Model 

The  flow  in  combustion  devices  is  known  to  be  predominantly  turbulent. 

To  account  for  this  turbulent  behavior  in  the  solution  of  the  time-averaged 
Navier-Stokes  equations,  a turbulence  model  is  introduced  to  define  an  effec- 
tive viscosity.  A review  of  turbulence  models  is  available  in  the  literature 
(see,  e.g.,  Refs.  19  and  20).  Prandtl  (Ref.  21)  was  perhaps  the  first  to 
introduce  a turbulence  model  when  he  postulated  that  the  time-averaged  shear 


If'.  Emmons,  H.  W.,  ed.:  Fundamentals  of  Gas  Dynamics.  High  Speed  Aerodyna- 

mics and  Jet  Propulsion,  Vol.  3 , Princeton  University  Press,  Princeton, 
N.  J.,  1058. 

11.  Launder,  B.  E.  and  D.  B.  Spalding:  Mathematical  Models  of  Turbulence. 

Academic  Press,  London,  1972. 

20.  Harlow,  F.  II.,  ed.:  Turbulence  Transport  Modeling.  AIAA  Selected  Re- 

print Series,  Vol.  XIV,  1973. 

21.  Prandtl,  L. : Bcricht  Uber  Untersuchingen  Zur  Aus geb i lde ten  Turbulenz. 

ZAMM.Vol.  5,  1925,  p.  136. 
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stress  and  the  time-avera.  ••  i / d.  -ity  -radie:.-  ax-  ; r j[  -rtionai  a.-  in  laminar 
fl aw,  and  that  the  engt  teal  the  - ■■  Length  whicl  enter:  th< 

relat.ion.'ii ! : Ls  proi  irti  nal  ( the  tur  it  r ton  thickness. 

Irandtl':;  mixing  iengtl  mod:  t been  ej  yed  iccessfully  by  a number  of 

investigators  (e.g.,  Refs.  It  • rariet;  of  pr  ■ emu  primarily  in- 

volving turbulent  flow  along  walls  an!  n free  tur: olent  flows.  A : i s advan- 
tage of  the  mixing  length  model  is  that  is  is  an  equilibrium  model  (:.e.,  • ur- 
bulence  is  assumed  to  be  produced  and  dissipated  locally)  and  it  requires  an 
ad  hoc  mixing  length  distribution.  Some  of  the  shortcomings  of  the  mixing 
length  model  have  been  overcome  for  many  cares  of  interest  by  the  introduction 
of  various  inultiequatlon  transport  models  of  turt  Hence.  A principal  disad- 
vantage  if  the  b L<  piatlon  • irl  ilenci  model  fr  . a nq  itat  i nal  viewpoint 
is  the  necessity  if  . ilvin  additional  trans]  srt  equations.  In  view  of  the 
preliminary  nature  of  the  present  w or  a mixing  length  model  is  employed  for 
several  reasons . First  >f  i . , ai  ad  : mlxi  . ength  :an  be  assumed  which 
does  give  a reasonable  resresentat 1 or:  of  • ne  t ,ri  .dent  process , ana  thus 
is  expected  to  give  n as  mal  Li  th<  retj  : re  li  ••  >ns.  S<  condly.  It  is  a 
definite  advanta.a  to  ke<  tl  :al  modi  s (t  irl  »lenc<  , chemistry,  drop- 
let, radiation)  relativi  y simple  to  a /.■.  rifi  sati  in  if  the  ba;  Lc  three- 

dimensional  numerical  calculation  procedure  at  this  point  in  time. 

die  formulation  of  the  mixing  length  turbulence  model  employed  In  thi. 
analysis  is  .,aced  on  the  mixing  length  distrit ution  suggested  by  Williamson 
(Ref.  25)  for  flow  in  ducts  and  nines.  It,  mathematical  terms  the  expression 
for  the  effective  viscosity  takes  the  form  (Ref.  27) 


^eff 

Re 


Pi 


2 ( 2e 


(12) 


where  the  mixing  Length  l is  given  by 


22.  Patankar,  S.  V.  and  D.  B.  Spalding:  Heat  and  Macs  Transfer  in  B mdar 

layers.  Intertext  Books,  London,  1 '70. 

27.  Maine,  G.  and  I!.  McDonald : Mixing  Lengtl  and  Kinematic  Eddy  Viscosit) 

a Compressible  Boundary  Layer.  AIAA  Journal,  Vol . • , 196R,  np.  /’s-hV. 
24.  McDonald,  H.  and  F.  J.  Camarata:  An  Extended  Mixing  Length  At  : reach  for 

Computing  the  Turbulent  Boundary  Layer  Development.  Proceedin'-.'  of  tin- 
AFOSR-IFi' '■-Stanford  Conference  on  Boundary  Layer  Prediction,  1 » . 

15.  Williamson,  J.  W. : An  Extension  of  Prandtl * s Mixinj  hen  tl  Che  ty. 

Applied  Mechanics  and  Fluids  Engineering  Conf-rence,  A;' ME,  Tune  1 » ■». 
Lilley,  p.  G. : Prediction  of  Inert  Turbulenl  Swirl  Flows.  AIAA  ■ v • r 

Wo.  72-699.  AIAA  5th  Fluid  and  Plasma  Dynamics  Conference,  June  l1!’. 
27.  Beer,  J.  M.  and  N.  A.  Chigier:  Combustion  Aerodynamics.  Wiiey, 

New  Voi-k,  1972. 


(13) 


77  = oi4(tH  e»p('-3r) 

Here  rQ  is  the  distance  from  the  wall  to  the  centerline  of  the  duct  or  pipe, 
and  y is  the  distance  from  the  point  in  question  to  the  nearest  wall. 

Special  consideration  must  be  given  to  calculation  of  turbulent  flow  in 
the  vicinity  of  walls  in  view  of  the  large  flow  gradients  which  occur  there. 
Since  the  expense  of  simply  increasing  the  number  of  grid  points  near  a wall 
may  be  considerable,  an  analytical  wall  function  formulation  has  been  employed 
in  the  present  study.  A set  of  three  universal  velocity  profiles  (Ref.  28) 
are  employed  in  the  wall  region,  corresponding  to  the  laminar  sublayer  (y+s4), 
a transition  region  (4  < y+  < 26),  and  the  logarithmic  law-of-the-wall  region 

(y+  * 26): 


where 


and 


♦ ♦ 
u = y 


u*  = c,ln  (l  + y*)  + c2 


u*  = Cj  Iny*  + c2 


for  y*  < 4 

+[(i-c,-c2  a)y*-  c2]  e~°y 
for  4 < y*  < 26 
for  y*>  26 


u = 


y* 


= y 


u* 


u 


# _ 


1/2 


(l4a) 


(l4b) 

(l4c) 


(15) 


(16) 


(17) 


28.  Walr. , A.:  Boundary  Layers  of  Flow  and  Temperature . The  M.I.T.  Ihress, 

Cambridge,  Massachusetts,  196 9. 
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In  Eq.  (15)  ~ denotes  the  total  velocity  component  parallel  to  the  vail,  and  u* 
is  the  nondimens ional  friction  velocity,  Eq.  (17).  The  finite  difference  form 
of  the  velocity  gradient  at  the  grid  point  adjacent  to  the  wall  point  's 
specified  consistent  with  the  appropriate  universal  profile  'iven  above.  In 
the  present  formulation  the  assumption  of  constant  shear  stress  the  immedi- 
ate vicinity  of  the  wall  is  utilized.  The  specification  of  the  velocity 
gradient  at  the  grid  point  adjacent  to  the  wall  (where  the  velocity  is  known) 
is  equivalent  to  imposing  a "slip"  velocity  at  the  wall  itself.  Hence  flow 
resolution  in  the  region  very  near  the  wall  is  sacrificed  to  attain  accuracy 
in  the  central  region  of  the  flow  field  where  the  combustion  processes  are 
concentrated.  However,  if  accurate  calculations  in  the  wall  region  are  re- 
quired at  a later  date,  refinements  to  the  wall  function  approach  may  be 
implemented  easily  in  the  present  computational  procedure. 

Chemistry  Model 

The  method  employed  for  introducing  the  effects  of  hydrocarbon  chemistry 
into  the  calculation  procedure  is  based  on  solution  of  the  fuel  conservation 
equation  with  a specified  pseudo-kinetic  chemical  rate  term.  The  pseudo- 
kinetic  approach  simplifies  the  computational  problem  of  incorporating  the  hydro- 
carbon chemical  energy  release  into  the  present  procedure,  since  local  chemical 
equilibrium  is  achieved  by  increasing  the  chemical  rate  constant  as  a function 
of  time  to  a sufficiently  large  value.  Furthermore,  a real  (global)  hydro- 
carbon kinetic  mechanism  may  be  incorporated  into  the  present  framework  in 
a conceptually  straightforward  manner . The  kinetic  nitric  oxide  (NO)  chemistry 
model  used  herein  assumes  that  NO  is  a trace  species  which  has  a negligible 
effect  on  mixture  properties  and  the  energy  deposition  in  the  flow  field. 
Therefore,  the  NO  species  conservation  equation  may  be  solved  separately  after 
determination  of  a steady  solution  with  hydrocarbon  chemistry. 

Nitric  Oxide  Chemistry  Analysis 

Solution  of  the  species  Eq.  (2)  to  account  for  convection,  diffusion, 
and  production  of  nitric  oxide  (lIO)  requires  an  expression  for  the  rate  of 
creation  of  nitric  oxide,  rpo,  due  to  chemical  reactions.  To  develop  this 
expression  knowledge  of  the  reaction  mechanism  by  which  nitric  oxide  is  formed 
is  required . A senerally  accepted  reaction  mechanism  for  NO  formation  and 
decomposition  in  post  flame  gases  is  that  proposed  by  Lavoie,  et  al.  (Kef.  2‘b'. 
It  consists  of  the  following  six  reactions: 


2'j.  Lavoie,  ;.  A.,  J.  B.  Heywood,  and  j.  C.  Keck:  Ejcjx  rimentai  u I 

Theoretical  Study  of  Nitric  Oxide  Formation  in  Tr.ternaJ  ’oirn'ustion 

Engines . Combustion  Science  and  Technology,  Vol.  1, 
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N2  + 0 cn  NO+N  (18) 

N + 02  = NO  +0  (19) 

N + OH  — NO  + H (20) 

H + N20  = N2  + 0H  (21) 

0 + N20  — N2  + 02  (22) 

0 + N20  = NO  + NO  (23) 


The  first  two  reactions,  Eqs.  (l8)  and  (19)  form  the  Zeldovich  mechanism 
(Ref.  30)  which  is  considered  to  be  the  principal  nitric  oxide  formation 
reaction  mechanism.  The  two  reactions  together  with  the  third  reaction,  Eq. 

(20)  which  assumes  minor  importance  under  fuel  rich  conditions,  form  the  ex- 
tended Zeldovich  mechanism  which  is  employed  in  the  present  study.  At  low 
temperatures,  when  nitric  oxide  concentrations  are  much  greater  than  equilibri- 
um values,  the  fourth,  fifth  and  sixth  reactions,  Eqs.  (21)  trhough  (23),  in- 
volving N20  as  an  intermediary,  may  become  important;  however,  because  overall 
reaction  rates  are  so  low,  the  net  effect  of  these  reactions  is  probably 
negligible.  Therefore,  these  reactions  have  been  neglected  in  the  present 
study.  The  results  of  Bowman  and  Seery  (Ref.  31),  as  well  as  other  investiga- 
tors (Ref.  32)  in  investigation  of  nitric  oxide  formation  kinetics  in  combus- 
tion processes,  lend  support  to  this  approach. 

With  the  reaction  mechanism  for  nitric  oxide  defined  by  the  extended 
Zeldovich  mechanism,  Eqs.  (18)  through  (20),  it  becomes  possible  to  develop 
an  expression  for  the  rate  term,  ruo>  entering  in  the  chemical  species  equa- 
tion, Eq.  (2).  The  mathematical  development  of  the  rate  term  is  based  on  the 
"Law  of  Mass  Action,"  which  states  that  the  rate  of  chemical  reaction  is  pro- 
portional to  the  active  masses  of  the  reacting  materials.  Thus,  referring  to 
Eqs.  (18)  through  (20),  rate  expressions  for  nitric  oxide  and  nitrogen  can  be 
written 

30.  Zeldovich,  Ya.  B.,  P.  Ya.  Sadounikov,  and  D.  A.  Frank-Kamenetskll: 

Oxidation  of  Nitrogen  in  Combustion.  Academy  of  Sciences  of  USSR,  Insti- 
tute of  Chemical  Physics,  Moscow- Leningrad,  19U7. 

31.  Bowman,  C.  T.  and  D.  J.  Seery:  Investigation  of  NO  Formation  Kinetics  in  the 

Combustion  Process:  The  Methane-Oxygen-Nitrogen  Reaction.  Emissions 

from  Continuous  Combustion  Systems.  Plenum  Publishing  Company,  New  York, 
1972. 

32.  Caxetto,  L.  3.,  L.  H.  Muzio,  R.  T.  Sawyer,  and  E.  S.  Starkman:  The  Role 

of  Kinetics  in  Engine  Emission  of  Nitric  Oxide.  Combustion  Science  and 
Technology,  Vol.  3,  1971. 
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k|fc  co  + k2,cNc  +k3fcNc( 


“klbl  NOrN  *'2bf'N0C0~  k3D  CNOCH 


= k|f  C N2C0  + k2t>LN0  C0+  **  3b  CNOCH 
'''lb  CN0  CN~*'2fCN  ^Og  ~ ^ 3f  CNCOH 

where  Ch  denotes  the  concentration  of  species  i (moles  per  unit  volume),  and 
k^.  and  k^  are  the  forward  and  backward  rate  constants  (volume/mole-sec),  res- 
pectively, for  reactions  (l8)  to  (20),  (subscripts  1,  2,  3).  Because  the 
relaxation  time  for  Eq.  (25)  is  several  orders  of  magnitude  shorter  than  for 
Eq.  (24),  it  is  a good  approximation  to  assume  a steady  concentration  for  CT^T. 
Setting  dCjj/dt  = 0,  Eq.  (25)  may  be  solved  for  C^: 


K|fCN2C0  + kg  b ^no  co  + k3b  CN0  CH 

^lbCNO  + ^2fC02+^3f  C OH 


This  equation  may  be  substituted  into  Eq.  (24)  to  yield  an  expression  for 
dCN0/dt  as  a function  of  the  concentrations  of  02,  N^,  0,  H,  OH,  and  NO,  and 
the  rate  constants  which  are  functions  only  of  temperature.  The  correct  non- 
dimensional  expression  for  the  rate  term,  r^Qj  to  be  employed  in  Eq.  (2)  is 


where  rjjo  is  the  mass  rate  of  production  of  nitric  oxide. 

The  rate  constants  governing  the  nitric  oxide  reaction  scheme,  Eqs.  (18) 
through  (20),  and  utilized  in  Eqs.  (24)  and  (25),  are  of  the  modified 
Arrhenius  form 


k.  = A.T  I exp  (-Bj/RT) 


where  data  for  the  activation  energies,  Bj,  the  frequency  factors,  Aj,  and 
exponent  Nj,  are  taken  from  Refs.  (33)  and  (34)  and  tabulated  in  Table  I.  The 
concentrations  appearing  in  Eq.  (24)  are  determined  using  the  chemical  equili- 
brium computational  analysis  of  Brinkley  (Refs.  35,  36),  with  the  amounts  of 
hydrocarbon  fuel  and  oxygen  available  for  reaction  determined  from  the  pseudo- 
kinetic  chemistry  calculation  described  below. 

It  should  be  noticed  that  time-mean  concentrations  and  temperature  are 
employed  in  Eqs.  (2k),  (26)  and  (28),  even  though  the  effects  of  turbulent 
fluctuations  in  these  quantities  may  lead  to  significant  errors  (Ref.  37)  in 
the  local  rate  of  production  of  nitric  oxide,  Eq.  (27).  The  implications  of 
these  effects  are  not  well  understood  at  the  present  time,  and  few  theoreti- 
cal modeling  approaches  to  the  problem  are  available  (e.g..  Ref.  SS'1  while 
even  fewer  have  been  evaluated. 

Pseudo-Kinetic  Hydrocarbon-Air  Chemistry  Analysis 

Tiie  pseudo-kinetic  chemistry  model  provides  a convenient  means  for 
introducing  the  effects  of  hydrocarbon  combustion  into  the  present  time- 
dependent  computational  procedure.  The  basic  chemistry  model  assumes  that  the 
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Combustion  Processes,  Section  C.  High  Speed  Aerodynamics  and  Jet  Propul- 
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the  combustion  process  ei  orgy  release  is  adequately  described  by  a single 
step  reaction  of  the  form 


Cf)  + ^2^2  ^ 


(29) 


In  the  reaction  (29)  the  hydrocarbon  fuel  CnHm  combines  with  Og  to  produce 
products  of  combustion  P,  with  the  nitrogen  in  air  being  treated  as  an  inert 
species.  The  necessary  species  conservation  equations  are 


d(Pm,) 

~a~P 

d{P< t>) 


_V-(pum,)+  — V + r, 


= — V -(pu<t>)  + J-  V (rmV4>) 
d\  Ke 


(30) 

(1  i) 


where  itiq  is  the  mass  fraction  of  unburned  fuel,  and  $ is  a mixture  fraction 
defined  as 


<1>  = Rs  m(  - m2 


(32) 


Here  m2  is  the  mass  fraction  of  oxidizer  (O2)  and  Rs  is  the  stoichiometric 
ratio,  i.e.,  the  ratio  of  oxidizer  mass  to  fuel  mass  for  reaction  (29): 

R = .W2  (33) 

S Wfuel 

where  W Oo  and  Wqueq  are  the  molecular  weights  of  oxidizer  (Og)  and  fuel  (CnH^), 
respectively,  and  v2  is  the  number  of  moles  of  oxidizer  in  reaction  (29). 

It  can  be  shown  that  the  mass  fraction  of  inert  species  (primarily  Np ) is 
given  by 

(3*0 


I-  f 


f0  + Rs 


— (Rs* 


where  f is  the  mass  fraction  of  O2  in  the  oxidizer  (fQ  = 0.2322  for  air^ 
mass  fraction  of  the  products  of  combustion  is  then 


The 


= i - m,  - m,  - m. 


(35) 


For  simplicity  it  has  been  assumed  that  the  reaction  (2°)  for  methane  proceeds 
.to'chiometrically,  so  that  vo=2  and  the  products  ..f  combustion  art  P^COp+’hv  ■ 
This  assumption  adequately  1 ;p resents  the  energetics  of  m< tl  ant -air  c ml  . ' 1 . 

• : resent  analysis  the  pseudo-chemical  production  rate  r-i  in  :■ 

:h  »si  r to  force  thi  fuel  mass  fraction  mq  to  its  equilibrium  va  le,  n . 

ri  = ^hc  (mi  “ mie) 
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The  rate  constant  kqr  is  increased  as  a function  of  time  to  a sufficiently 
large  value  to  insure  achievement  of  equilibrium.  The  justification  for  the 
chemical-  equilibrium  approximation  is  based  on  the  observation  that,  for 
many  combustors,  the  temperature  and  pressure  conditions  are  such  that  the 
fuel  oxidation  reactions  go  to  completion  rapidly.  A nonequilibrium  chemis- 
try model  could  be  implemented  in  the  present  procedure  at  a later  time  if  it 
is  warranted  by  comparison  with  experimental  data. 

For  present  purposes  it  is  adequate  to  assume  that  at  a given  point  either 
all  the  fuel  or  all  the  oxidizer  is  consumed  in  the  equilibrium  state.  Hence, 
for  combustion  governed  by  reaction  (29)  the  equilibrium  fuel  mass  fraction  :s 
given  by 


miP  = 7T-  & 

‘ ® R c 

0 

Al 

■e 

0 

1 ' S 

(37) 

O 

11 

a> 

E 

for  <t>  < 0 

The  energy  deposition  rate  term,  £ vAx . , in  the  stagnation  enthalpy  equation 
(5)  takes  the  following  form  for  1the  combustion  reaction  (2n): 

Zr,h!:  rt  [h|-hp  +Rs(h2'hp)]  ( 

i 

f f f 

where  h-j_,  h2,  and  h are  the  heats  of  formation  for  the  fuel,  oxidizer  and 
combustion  products , respectively. 

The  species  concentrations  which  are  required  for  the  NO  chemistry,  Eq. 
(24),  are  determined  using  a chemical  equilibrium  analysis  (Ref.  35,  -'<>). 
However,  the  present  application  may  be  considered  as  a "parti ad- equilibrium" 
analysis,  since  the  amounts  of  hydrocarbon  fuel  and  oxygen  which  have  reacted 
are  determined  from  the  pseudo-kinetic  model  as  follows.  The  solution  to 
Eqs . (30)  and  (31)  yields  the  mass  fractions  of  (uribumed)  fuel  and  xygei  . 
m^_  and  m,,  existing  at  a given  point,  [t  easily  shown  that  the  • ta.  fuel 
and  air  mass  fractions  present  in  the  air  Mice  if  cl  ejnj  :al  r acl  i ni  r ] 
by 


m 


1 


$ +f0 

Rs  + f 0 


(3") 


and 


m 


]ir 


r I - m 


-»  + Rs 
Rs+f0 


(40) 


It. 


fraction  is 


••>  • it  the  corresi  ndinr  oxygen  mass 

m2  r f o mau  (4l) 

Finally,  t’nc  mis  fraeti  n .f  fuel  and.  oxy  en  (ioq^)  which  have  reacted 

are 


mHC  = m|—  m, 


(42^ 


m 


V 


m ^ — m p 


(43) 


m . • . which  disai  peared  kinetically 

Lbriui  manner.  In  ‘-he  partial  equilibrium  calculation 
..  ri  from  • cies  not  available  for  reaction,  m^  and  mo,  are 

inert  wil  equi  riui  ihemist  rji  r ledure.  [heir  1 

:•  juilibrium  composition  come  mly  from  their  contribution 

ter  j ■ ••  . is<  of  the  above  formulation  with  the  con- 

straints f specified  ■ res sure  and  enthalpy  permits  the  calculation  of  the 

ippearii  n i«  (24) . 


. '-•);lf  : ' ' i.  . 


' • eel  f liqid  s'  jombustor  mizi  rs  results  1 n 

•.he  forma: ion  f a ■ at : a ' die  - r'sution  of  liquid  droplets  within  a specified 
range.  ■ *ation,  ration,  and  < / ntuaJ  I un lin  of  these  in >;  - 

m tant  I rf  rmance  tf  lombustors . According- 

. •>  , U" ; ■ .)•  hayed  in  tlr's  study  to  account  for  tiie  ; nfluence  of 

dro;  s.'s  n • lEii’U  r : -rf<  nuance. 


il.  ••out  ; in  ir  >:*.let  spray  is  lirectly  influenced  by  the  droplet 
re  s'  , • . ■ T.  '5  and  125  microns  for  ::ust  practical  liquid 

iditi  lr  size,  dro:  Let  combustion  rates  are  al. 

•i  i’1'  i i 1 "s'  !■  : <•  .[•<•  us  inert ' and  the  relative  velocity  between  the 

: exity,  the  e ff<  its  »f  these 

imin  *at  avi  nly  beej  Li  \n  s to  ited  for  sin 
ited  dr  . r . t ibsequent  iroplet  model  develo]  mi  at  Is 

■ ■ eocti itributiom  :'r  i sin  L<  droj  el 

analyses . 

: s;  ic count  for  c mn  sett  1 . 

. . *a1  imi  Lquid  droplets.  Che  model  is  based 
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solution  of  the  particle  fraction  equations  (Refs.  39  and  40)  having  as 
dependent  variable  the  particle  fraction,  f^,  which  represents  the  mars 
fraction  of  liquid  contained  in  an  incremental  droplet  size  range  (see  Fie. 

2).  Tiie  particle  fraction  equation  (3)  is  again 

= -wof,)  * (r0vi,|.R 

where  Tp  is  the  liquid  phase  eddy  diffusivity  coefficient,  and  K-j  represents 
the  rate  of  production  of  particles  in  class  j due  to  all  sources.  One  par- 
ticle fraction  equation  is  needed  for  each  particle  class  studied.  Although 
in  principle  the  method  allows  for  an  infinite  number  of  particle  classes, 
computer  time  and  storage  considerations  have  set  a limit  of  four  particle 
classes  in  the  present  analysis.  Solution  of  the  particle  fraction  equations 
is  sufficient  to  define  the  distribution  of  droplet  sizes  existing  at  each 
location  in  the  flow  field. 

Implicit  in  the  use  of  the  droplet  model  are  certain  underlying, 
assumptions  which  are  listed  as  follows: 

(1)  The  droplets  form  a cloud  of  suspended  particles. 

(2)  The  volume  of  droplets  is  negligible  compared  to  gas  volume. 

(3)  Relative  velocity  of  droplets  and  surrounding  gas  is  negligible. 

(4)  Temperature  and  density  of  droplets  is  uniform. 

(5)  Gas  and  liquid  phase  transport  coefficients  are  equal. 

As  a consequence  of  assumption  (1)  droplet-droplet  interaction  is  neglected. 
Assumption  (2)  is  merely  a statement  that  droplet  density  is  much  greater  than 
gaseous  density,  which  simplifies  the  governing  equations  and  associated  c im- 
putations. Assumptions  (3)  and  (4)  also  lean  to  a simplification  >:'  the 
governing,  equations  because  under  these  conditions  only  terms  describing 


39.  Jibs  M.  M.  and  B.  !•'.  Morgan:  Mathematical  Model  of  Combustion  o' 

Solid  Particles  in  a Turbulent  Stream  with  Recirculation.  J rarn&l 

the  Institute  of  Fuel,  December  1°70. 

40.  Spalding,  D.  B. : Mathematical  Models  of  Continuous  Combustion.  Emission 

from  Continuous  Combustion  Systems,  Plenum  Publishing  Company,  New  York, 
1972. 
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interphase  mass  transfer  need  be  considered.  Assumption  (5)  eliminates  the 
need  for  defining  liquid  phase  transport  coefficients  about  which  little  is 
known.  Justification  for  assumptions  (3)  through  (5)  come  from  '.he  fact 
that  droplet  sizes  in  typical  combustors  can  generally  be  expected  tc  be 
small  both  in  regard  to  physical  size  (low  micron  range)  and  witn  respect  to 
the  microscale  of  the  turbulence.  Under  these  conditions,  the  droplets  re- 
spond closely  to  the  mean  and  fluctuating  components  of  the  gaseous  motion 
and  can,  therefore,  be  expected  to  exhibit  similar  transport  behavior.  (Prop- 
let  sizes  large  in  comparison  to  the  turbulence  microscale  do  not  follow  the 
gas  motion  but  behave  as  though  suspended  in  a laminar  flow  field  having  the 
same  mean  motion  as  that  of  the  turbulent  flow.)  Under  these  conditions 
there  is  not  droplet  diffusion.  Solution  of  the  particle  fraction  equation, 
Eq.  (3)j  is  straight forward  in  the  computational  procedure  subject  to  appro- 
priate expressions  for  the  boundary  conditions  and  the  source  term.  The 
boundary  conditions  are  readily  formulated  for  a typical  combustor  (see  sec- 
tion on  Boundary  Conditions).  The  source  term  in  Eq.  (3)  is  evaluated  for 
particles  in  class  (j)  by  assuming  that  mass  gain  is  due  to  vaporization  of 
particles  originating  in  class  (j+1)  which  thereby  enter  class  (j).  Mass 
loss  for  class  (j)  is  due  to  vaporization  of  particles  in  class  (j)  which 
thereby  enter  class  (j-1).  That  is,  as  indicated  in  Fig.  2,  an  increase  in 
the  number  of  particles  within  a given  class  is  due  only  to  the  decrease  in 
radii  of  particles  which  are  initially  in  classes  having,  larger  radii.  Under 
this  study  an  increase  in  particle  fraction  due  to  condensation  from  the 
gaseous  phase  is  excluded  from  consideration.  By  treating  the  loss  or  gain 
in  particle  fraction  due  to  the  change  of  radii  from  one  particle  class  to 
another  as  a movement  of  particles  whose  radii  are  changing  at  a rate  dr/dt, 
the  rate  of  change  of  particle  fraction,  f j , at  the  upper  and  lower  radii 
boundaries  can  be  formulated  as  follows: 


Rate  of  movement  of  particles 
from  above  rjH  ^ to  below  r j q 


Pi  r' /dr\ 

r _ r \ dWi*1 

j+2  )+l 


Rate  of  movement  of  particles 

from  above  r . to  below  r ; 

J 


(WO 


5) 


In  addition  to  the  two  terms  represented  by  Eq:'- . (410  and  (45),  which  repre- 
sent mass  transfer  across  adjacent  particle  class  boundaries,  the  source  term 
in  the  particle  fraction  equations  must  contain  a third  term  which  represents 
loss  of  mass  directly  to  the  gaseous  phase.  This  term  can  lie  represented  as 


Loss  of  mass  to  gaseous  phase  f — - dr 

J r - r r dt 

l + i 


(4t  ^ 


J 


Employing  Eqs.  (44)  through  (46),  the  complete  expression  for  the  source  term 
for  the  j^h  particle  class  becomes 
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Assuming  that  the  droplets  are  spherical 
easily  be  shown  that 


and  of'  constant  density,  it  can 


dr  _ _ m 
dt  47rr2/°d 


(48) 


where  is  the  droplet  density  and  m represents  the  rate  of  mass  loss  of  a 
particle  due  to  vaporization.  A semi-empirical  expression  for  m 'used  in  t!  's 
study  has  been  derived  from  a single  droplet  analysis  (Ref.  4l  and  42).  and 
has  the  form 


m = Avr  md 


(4d) 


and 


Me  ff 


md  = 


Pr 


eff 


In  {i  + hv[cp(Tg-Tb)+  hc^i3-]| 


(50) 


w'here  hy  is  the  heat  of  vaporization  per  unit  mass  for  liquid  fuel  at  the 
temperature  T (the  surrounding  gas  temperature),  Tb  is  the  liquid  boilinv 
temperature,  hc  is  the  heat  of  combustion  per  unit  mass  of  oxygen  for  the 
chemical  reaction  considered,  s is  ratio  of  oxygen  mass  to  fuel  mass  for  the 
reaction  considered,  and  im,  r is  the  mass  fraction  of  oxygen  in  the  surroundin' 
gas. 


Using  Eq.  (49)  in  Eq.  (48),  one  obtains 


dp  _ _ m d 

dt  ’ Pd  r 


(5D 
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Physical  Parameters  on  the  Burning  Rate  of  a Liquid  Droplet.  Fifth 
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Substituting  this  expression  into  Eq.  (4-7)  and  integrating  the  last  term 
of  that  equation  yields  the  final  form  for  the  droplet  source  term: 

RPr_^md  (_{)_[ i +J 111! ) (5 

1 '°d  | ^ S+l  ^ Vl  fJ+l(r|  + 2 Vl)  J 

Equation  (52)  represents  the  rate  of  fuel  loss  to  the  gaseous  phase  as  a 
result  of  vaporisation,  it  also  includes  the  mass  burning  rat*  f the  trip- 
let because  of  the  equilibrium  assumption  employed  in  the  analysis  wherein 
the  chemical  reaction  is  presumed  to  go  to  completion  and  to  occ^ir  instanta- 
neous ly. 

Radiation  Model 

It  is  well  known  that  gases  and  other  substances  at  high  temperatures 
emit  energy  In  the  form  of  electromagnetic  radiation.  The  purpose  of  the 
radiation  model  employed  under  this  investigation  is  to  define  the  radiant 
energy  contribution  to  combustion  heat  transfer  rates.  Specifically,  this 
model  is  employed  to  define  the  radiant  enemy  source  term,  -V • q^,  entering 
the  energy  equation  (5),  through  which  coupling  between  the  radiation,  con- 
vection, and  conduction  modes  of  heat  transfer  occurs.  Unless  the  gaseous 
medium  in  a combustion  chamber  is  strongly  absorbing,  the  effect  of  radiation 
will  -enerally  be  to  reduce  the  peak  temperatures  in  the  flow  field  while  wall 
surface  temperatures  will  increase  due  to  the  possibly  large  radiant  heat 
flux  reaching  the  surface.  In  these  cases,  accurate  prediction  of  combustor 
performance  and  emission  characteristics  will  require  a reasonable  representa- 
tion of  the  radiative  transfer  process.  The  radiative  energy  transfer  pro- 
cess is  inherently  different  from  conductive  and  convective  heat  transfer 
processes  since,  in  the  latter  two  cases,  enemy  is  transferred  by  molecular 
collision  and  transport,  whereas  radiative  enemy  transfer  does  not  require 
molecular  contact.  The  emission  and  absorption  of  thermal  radiation  is  due  to 
■ transitions  between  the  energy  levels  of  the  atoms  or  molecules  of  the  gas, 
and  transitions  involving  free  electrons. 

The  preliminary  radiation  transport  model  employed  in  this  study  is  a 
discrete-flux  model  (Ref.  il).  Similar  methods  have  been  employed  with  some 
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success  by  other  investigators  (e.g.,  Refs.  44  through  46),  and  they  are 
described  in  the  literature  (Refs.  46  through  48).  The  procedure  developed 
in  Ref.  43  is  an  ad  hoc  extension  of  the  one-dimensional  transport  model: 
proposed  by  Schuster,  Schwarzchild,  and  Hamaker  (see  Ref.  48  for  other  refer- 
ences), in  which  discrete  radiation  fluxes  in  each  of  the  positive  and  nega- 
tive coordinate  directions  are  considered.  In  the  present  study  the  radiation 
field  will  be  assumed  to  be  nearly  axisymmetric  for  axisymmetric  combustor 
geometries  so  that  only  fluxes  in  the  radial  and  axial  directions  need  be  con- 
sidered (the  four- flux  model,  Ref.  43).  For  rectangular  geometries  a six-flux 
model  (Ref.  14)  is  easily  implemented.  Since  the  four-flux  model  has  been 
described  elsewhere  (Refs.  43  and  49),  only  the  six-flux  model  will  be  out- 
lined here. 

t 

Consider  forward  and  backward  fluxes,  qq  and  q^ , in  the  coordinate 
direction  xq.  Hie  governing  equations  for  a grey  emitting  and  absorbing  me- 
dium in  local  thermodynamic  equilibrium  including  isotropic  scattering  may  be 
written  in  Cartesian  coordinates  (i  - 1,  2,  3)  as: 


<“o  ♦«.><:  + (-^)<v4*  f-z  (»[♦«;) 


■$*-=-  {%*%)  i?  ♦ 


(?)  f ! K-V> 


where  o'-  and  og  are  the  so-called  "flux"  absorption  and  scattering  coefficients 
(Ref.  43),  a is  the  Stefan-Boltzmann  constant,  and  T is  the  absolute  tempera- 
ture. Before  solution  of  the  transport  equations  (53)  and  (54)  can  be 
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attempted  the  mean  absorption  coefficient  cv^  and  scattering  coefficient  a„ 
mast  be  specified.  Because  scattering  is  generally  of  secondary  importance 
in  a particle- free  gaseous  medium  of  the  type  under  consideration  in  the 
present  study,  the  scattering  coefficients  will,  not  be  defined  in  this  analy- 
sis. For  an  evaluation  of  the  absorption  coefficients,  reliance  is  placed 
on  the  data  available  on  emissivity  of  gases  and  gas  mixtures  (Refs.  50  to 
52).  The  true  absorption  coefficients  for  air,  water  vapor,  and  carbon  diox- 
ide have  been  presented  by  Cess  (Ref.  53)  using  the  emissivity  data  available 
from  the  literature.  For  the  present  purposes  it  will  be  assumed  that  the 
necessary  (frequency  averaged)  absorption  coefficients  are  known  functions  of 
the  thermodynamic  state  (e.g.,  pressure  and  temnerature) , and  will  not  be 
considered  further. 


A considerable  simplification  of  the  transport  equations  (53)  and  (54) 
may  be  realized  Dy  writing  them  as  a system  of  sc-ccnd-order  equations,  so 
that  only  three  additional  equations  must  be  solved  rather  than  six.  For  this 
purpose  the  following  definitions  are  needed: 

Gi  = iK  + ^)  (55) 

Qj  = q*  - q[  (56) 

Manipulating  Eqs.  (53)  and  (54)  with  Eqs.  (55)  and  (56)  yields  (with  the 
scattering  coefficient  retained  for  generality): 
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2 dG, 
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results  of  both  the  Schuster-Schwarz child  approximation  (isotropic  u.-ular 
distribution  of  intensity)  and  the  Schuster-Hamaker  approximation  (unidirec- 
tional radiation ^ reveals  that  the  flux  absorption  coefficient  (cO  should 
equal  da,  and  a^  in  these  two  cases,  respectively  (where  aa  is  the  true  ab- 
sorption coefficient).  Since  information  regarding  the  angular  distribution 
of  intensity  is  not  known,  it  is  necessary  to  consider  the  sensitivity  of  the 
predictions  using  this  model  both  when  performing  calculations  and  when 
assessing  the  accuracy  of  the  model. 


Boundary  Conditions 

Solution  of  the  governing  partial  differential  equations  requires 
specification  of  boundary  conditions  for  each  of  the  equations  on  all  boun- 
daries of  the  computational  region.  At  a solid  wall  the  density  is  determined 
.implicitly  using  a three-point  one-sided  difference  approximation  of  the  con- 
tinuity equation.  The  tangential  velocities  at  a wall  are  determined  impli- 
citly through  the  use  of  the  wall  function  approximation.  The  boundary  con- 
ditions for  each  equation  are  given  below  in  terms  of  the  normal  and  tangen- 
tial coordinates  (n,  s^,  S2)  at  each  surface.  The  radiation  boundary  condi- 
tions are  cons  ide  red  separately  at  the  end  of  this  section. 
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■ exil  i am  boundary  conditions  attempt  to  simulal  the  nearly  in 
tional  nature  of  the  flow  presumed  to  occur  there;. 
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where  un  and  arc  the  i.cmd.  radial  he;:  flux  ai;..;  ji.e  r .-lir-ct-  i 
ard  •;  . respectively.  Using  Eqs.  ( 

yields 
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where  n if  the  nona&l  coor  ii.nte  direct  ed  ciws;,'  t v m * ’•  •.*  .'uri’^.ce  1 y.  K .-\i  • • 
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In  this  case  there  is  no  emir-,  ion  or  reflection  so  that  ohc  incorl::. 
radiation  must  be  sped 'Med, 

qn  = «,n 

and  one  obtains 
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Tntr  duct ion 

Ex!  • : 1 ime-avi raged  Navier-Stokes  equations 

are  rare  due  to  their  high  order  arm  coupled  nonlinearity,  As  a result 
tm<  rici  thods  i Ly  be  employed  for  the  solution  of  these  equa- 
te ns.  in  addition,  one  of  the  ma  j >•  bsti  1<  he  lutioi  f th  ■ ree- 

i na]  >kes  equations  is  the  large  amount  of 

'input er  rime  required,  and  consequently  efficient  ccoputat ional  methods  are 
high] ;/  d s i rat  I /ious  method  foi  ample.  Refs.  5 to  • have 

• <;■>  based  on  explicit  difference  schom-  ■-  for  the  u - ;v  form  of  the  g vern- 

bjec  n - ■:  tabiJ  j • i i - 1 n > thi 

size  of  the  time  step  relative  to  the  spatial  mesh  size.  These  stability 
limits  usually  correspond  to  the  well-  • • • rid  - y : - 

itt-ior  and  in  some  methods  c...  an  nhiit.  nd  \\<~c  u.  stability  condition.  A 
hey  disadv  ntage  of  such  conditior.ally  stable  methods  is  that  the  maximum 
time  step  is  fixed  by  the  -:uatia]  mesh  siz  rathe:  than  th-  physical  time 
dependence.  If  a steady  soluti  n is  being  computed  as  tire  asymptotic  limit 
r large  timi  f ■ insti  y Lutioi  the  resei  : . :ti  . then 

using  a small  time  step  requires  a large  number  of  steps  to  reach  the  steady 
so] ut ion . 
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In  contrast  to  most  la. licit  methods,  many  implicit  methods  tend  to  be 
stable  for  large  tic.-  : tr,  . , ind  nence,  offer  the  prospect  of  substantial 
increases  in  -mputarional  efficiency,  provided  of  course  that  the  computa- 
tional effort  per  time  ste;  is  competitive  with  that  of  the  conditionally 
stable  methods.  An  implicit  method  based  • n an  alternating-direction-implicit 
(ADI)  scheme  has  been  developed  at  >iTRC  by  Briley  and  McDonald  (Ref.  15)  for 
solution  of  the  three-dimensional,  compressible  Navier-Stokes  equations.  This 
method  will  subsequently  be  referred  to  as  the  MINT  (i.e.,  Multidimensional 
Implicit  Nonlinear  Time -Dependent ) procedure.  The  MINT  computer  program  has 
provided  the  basic  numerical  framework  for  development  of  the  present  combus- 
tor flow  analysis.  The  MINT  method  has  been  described  in  detail  in  Ref.  15, 
but  it  will  be  summarized  here  for  completeness.  The  method  can  be  outlined 
briefly  as  follows;  The  governing  equations  are  replaced  by  an  implicit  time 
difference  approximation.  Terms  involving  nonlinearities  at  the  implicit  time 
level  are  linearized  by  Taylor  expansion  about  the  solution  at  the  known  time 
level,  and  spatial  difference  approximations  are  introduced.  The  result  is  a 
system  of  multidimensional  coupled  (linear)  difference  equations  for  the  de- 
pendent variables  at  the  implicit  time  level.  To  solve  these  difference  equa- 
tions, tiie  Douglas-Gunn  (Ref.  60)  procedure  for  generating  alternating-direc- 
tion-implicit  (ADI)  schemes  as  perturbations  of  fundamental  implicit  differ- 
ence schemes  is  introduced.  This  technique  leads  to  systems  of  one-dimen- 
sional coupled  linear  difference  equations  which  can  be  solved  efficiently  by 
standard  block-elimination  methods.  No  iteration  is  required  to  compute  the 
solution  for  a single  time  step,  and  for  a k-point  one-dimensional  difference- 
operator  only  2k-l  arrays  of  information  need  be  in  coie  at  any  one  time. 

The  scheme  is  therefore  very  economical  in  i erms  of  computer  memory 
requirements . 

difference  Notation 


The  governing  equations  (-5'  have  three  characteristics  which  are  given 
. pecial  consideration  in  the  numerical  formulation;  the  equations  are  (2  non- 
linear, ( ')  coupled,  and  (3  multidimensional . Before  proceeding,  some 
lifference  notation  is  in: i nlnoe-l . The  flow  egion  is  discretized  by  grid 
points  havirg  equal  spacing:: . Ax,,  AxjS  and  Ax.  in  the  Xj  , x and  x.  direc- 
tions, respectively,  and  an  arbi U ary" time  step.  At.  provisions  for  non- 
uniform  grid  spacing  will  be  introduced  subsequently.  The  subscripts  i.  . k 
and  superscript  n are  grid  p int  indices  an-  iciated  with  x-,  . x . and  x^.  and 
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t,  respectively.  Thus,  jiV  , k denotes  </>  (x.  , x x , t°)  where  ^ can 

represent  any  of  the  dependent  variables.  The  subscripts  are  frequently 

omitted  if  clarity  is  (reserved,  so  that  /“  is  equivalent  to  p . . For 

i , J • K 

convenience . the  following  shorthand  difference-  operator  notation  is  used  fr 
derivative  difference  formulas; 
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with  analogous  definitions  for  ,,  6",  6^  and  6'^.  It  is  assumed  that  the 
solution  is  known  at  the  n level,  tn,  and  is  desired  at  the  (n+1)  level.  t'‘  \ 


Linearization  Technique 

The  large  time-step  capabilities  of  implicit  methods  can  place  great 
demands  on  the  linearisation  technique  employed.  Indeed,  the  favorable 
stability  properties  of  implicit  methods  can  be  severely  compromised  by  an 
inadequate  linearization.  The  present  technique,  which  can  be  used  only  in 
conjunction  with  an  implicit  difference  scheme,  permits  the  implicit  solution 
of  coupled  nonlinear  equations  in  one  space  dimension  by  a one-step  noniter- 
ative scheme.  This  feature  is  retained  for  multidimensional  problems  by  using 
ADI  techniques.  The  linearization  is  accurate  when  variables  change  by  rela- 
tively small  ^mounts  during  a time  step,  and  consequently,  the  accuracy  of  the 
linearization  can  be  controlled  by  varying  the  time  step.  The  linearization 
technique  is  alsc  convenient  for  the  implicit  treatment  of  coupled  nonlinear 
boundary  conditions,  and  this  latter  feature  has  been  found  to  nave  a highly 
favorable  effect  in  the  stability  of  the  overall  method  (Ref.  15). 

A technique  is  described  for  deriving  linear  implicit  difference  approxi- 
mate ns  for  nonlinear  differential  equations.  Tne  technique  is  presented  with 
reference  t..  the  f .11  wing  first-  rder  equat.i  n in  one  dependent  variable, 

(x,  t): 
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Log  of  l ted  torn  is  as  yet  unspecified. 
Making  use  of  Chain-rule  differentiation,  the  bracketed  term  in  Eq.  (77)  is 

' ; * - ; th(  IfJ  - ■ ■ liffer- 

ences  and  centered  spatial  difference  t obtain  the  following  implicit  diier- 
ence  scheme: 


± n*i  , n 
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At 
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on  examination,  it  can  be  seen  that  Eq.  ( ")  is  linear  in  ' and  that 
all  other  quantities  are  either  known  or  evaluated  at  the  n level.  Because  of 
the  spatial  difference  operator,  6y,  Eq.  (78)  contains  , nd  /■'  ] : 

consequently,  the  system  of  linear  equations  generated  by  writ  . . 

a ■ 

he  implicit  ; n be  written  in  t 

' ' • cai  ■ . ■ ed  • ■ f 

niques  for  i riding. nal  systems  (see,  e.g. . Ref.  H ).  The  tridiagonal  matrix 
' r " • ’ ' i Bq.  in  the  i 
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When  a i : ! i - : r 


The  linearization  technique  will  to*  applied  t the  couple  i jjcvijr- 
lation  n th  fo  . . ttion,  but  Lrsl  • ei  bsei 

>priate.  liffet  theme  1 t-or  • y i time , 1 

bee  . constructed  fro-,  the  backward  difference  ret)-'  icement.  Eq.  (77 ». 

• a •••••!  •.  ra]  teuracy  can  be  increased  either  by  us  lug  a Crank-Nicolson 
terii  bout  tn  ...  f averaginj 

b“t.ver-n  the  n and  (n*-l  levels ) or  by  using  addition-.1  time  levels  ir  • 
difference  replacement  uf  time  del  ivative-?.  .1  w^ver,  t --rtporul  -.c.-i  • ••cy  i 
: t considered  iimr.  > i.ant  for  computing  toady  solatj  a.  a-  i : * a time-  .ft 
; therefore , t!  y backward  fon 

t i ■ n ■ : ■ ■ . . rning  equal 

) ha  conservation  fur.,.  (-■  .1  ) , then  the  i ine-ari  t . : } r.  -e.  prctvrv  . 
the  conservation  torn  for  steady  solutions. 


Application  of  the  Method 


The  extension  of  the  numerical  metho  : 1 . U.t  *.'<•  hv  nii.  n;  is  - red  : *v 
=■  1 1 *.-rtiu  ting-direction  implicit  (ADI  • technic-;.  . Th-  original  I moth  i wi- 
lt. r iuoed  by  1 eaceman  and  Rachferd  (f.e-f . e anti  i juntas  (Ref.  • _•  ; hew-- vt  r, 
1 • liter  ting-di  ction  ncept  ha  sine  beei  xpat  and  gen 
discussion  of  various  alternating-direction  technique*;  ir  given  by  Ki'eh-'I1 
(Re!  . 64  ••  nd  Yanenko  (Ref.  65) . Although  much  .•»  t.he  attention  1--.  *■  1 * 

ADI  method'  has  focused  on  single  line  r equal i • :.e*l  id  for  r.  nlineat 

ysten  : ggesl 
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The  present  technique  is  an  application  of  the  very  general  procedure 
developed  by  Douglas  and  Gunn  (Ref.  60)  for  generating  ADI  schemes  as  per- 
turbations of  fundamental  implicit  difference  schemes  such  as  the  backward- 
difference  or  Crank-Nicolson  schemes.  As  with  all  ADI  schemes,  Douglas-Gunn 
schemes  compute  the  solution  for  a time  step  in  a sequence  of  intermediate 
steps,  each  of  which  usually  involves  treating  derivatives  implicitly  in  one 
of  the  coordinate  directions.  Although  the  solution  of  multidimensional 
implicit  difference  equations  is  normally  time  consuming,  the  one-dimensional 
difference  equations  appearing  in  these  intermediate  steps  can  be  solved 
efficiently,  and  this  is  the  attraction  of  ADI  methods  in  general  over  other 
implicit  methods.  A less  obvious  feature  of  schemes  generated  by  the 
Douglas-Gunn  procedure  is  that,  since  rhey  give  the  solution  to  a fundamental 
implicit  difference  scheme  except  for  a small  perturbation  error,  they  can  be 
viewed  as  an  approximate  technique  for  solving  the  fundamental  difference 
scheme,  rather  than  as  difference  schemes  themselves.  This  unconventional 
view  of  the  Douglas-Gunn  procedure  provides  guidance  in  formulating  the  pre- 
sent method  by  indicating  that  the  governing  equations  can  be  linearized 
and  differenced  in  accordance  with  the  backward  difference  scheme  (or  another 
implicit  scheme)  before  the  alternating-direction  aspects  of  the  method  are 
introduced.  Douglas-Gunn  schemes  also  appear  to  have  an  advantage  over 
locally  one-dimensional  (LOD)  or  "splitting"  schemes,  whose  intermediate 
steps  do  not  satisfy  the  consistency  condition.  The  lack  of  consistency  in 
the  intermediate  steps  can  present  difficulties  in  the  treatment  of  some 
boundary  conditions  and,  according  to  Yanenko  (Ref.  65,  p.  33)  does  not  permit 
the  use  of  asymptotically  large  time  steps. 


The  numerical  method  is  essentially  an  application  of  the  linearization 
technique  of  the  previous  section  to  the  coupled  system  of  governing  equations 
(1-5).  These  equations  are  written  in  backward  time  difference  fora,  and  non- 
linear implicit  terms  are  linearized  by  expansion  about  tn.  The  viscous  force 


terms  in  Eq.  ("•)  which  contain  mixed  derivatives  (i.e.,  & 


dy  . for  i i j ) 


are  most  easily  treated  explicitly  by  evaluation  at  time  level  n.  Although 
mixed  derivatives  can  be  differenced  implicitly  within  the  Douglas-Gunn  frame- 
work, this  woiild  increase  the  number  of  intermediate  steps  and  thereby  com- 
plicate  th<  . >luti  n procedure.  For  the  solutions  presented  her<  and  test 
cases  computed  by  Briley  and  McDonald  (Ref.  15  , the  explicit  treatment  of  he 
aforementioned  viscous  terms  had  no  observable  adverse  affect  on  stability. 

optional  artificial  viscosity  terms  having  the  form 
) ^ n+l  are  added  to  the  difference  equations,  where 

respecti vel.v , in  the  con- 


In  three  dimensions 
/ n t n n 

(£!  6i  f e3 

/ represents  p,  u,  v,  w.  H,  mi, 

x - momentum,  energy,  fuel  species,  hybrid  species,  par* i 


tinuity,  x^-,  x -, 


f . . and  m„ 
J 


cle  fraction,  and  N1  . peci-  quat.ions.  The  difference  equations  obtained 
by  the  procedure  just  outlined  represent  '<  linearized  backward  difference 


V* 


scheme . The  equations  can  be  arranged  according  to  time  and  space  derivatives* 
and  written  in  the  following  matrix  operator  form  (Ref.  15): 


= n 
A 


' At  ' 


- n - n 

= D.  4> 


+ 


= n — n*i 
D2  <*>  + 


= n 7(i*i  _ n 

D.  + Sn 


(8o) 


-n 

Here  A is  a (mxm)  matrix  containing  the  time  derivative  coefficients,  where 
m is  the  number  of  equations  being  solved;  f is  the  column  vector  of  the 
dependent  variables;  D^  > D .. n , and  D-^n  are  (mxm)  matrices  containing  three- 
point  difference  operators  associated  with  the  coordinate  directions  x.,  x,, 
and  , respectively;  and  S is  a column  vector  containing  only  n-level 
terms . 

Since  the  multidimensional  implicit  system  with  coefficients  generated 
by  Eq.  (80)  is  difficult  to  solve,  the  Douglas-Gunn  (Ref.  60)  technique  is 
applied  to  Eq.  (80)  to  generate  an  ADI  scheme.  With  the  observation  that 
the  Douglas-Gunn  procedure  is  being  applied  to  a coupled  system  of  equations, 
the  following  three-step  scheme  is  obtained. 


= n_7*  =n,n  r-n  -rt 

D,  9 + D2  9 + D3  9 + S 


r(V-) 

»n(— fl|^  )-  W f + B;  i +B,+"-.sn 


= n I =n  r**  =n  7 1 
D,  9 + D?  9 + D, 


_n 
+ S 


(811 


(82) 


(83) 


where  <f>  , p , and  <f>  are  the  intermediate  solutions.  Note  that  at  each 
step  of  the  scheme,  one  more  coordinate  direction  is  treated  implicitly,  and 
that  the  most  recent  approximation  to  ji5  is  not  always  used,  as  this  would 
adversely  affect  the  stability.  Douglas  and  Gunn  were  able  to  show  under 
fairly  general  assumptions  that  the  .ADI  scheme,  Eqs.  (81-83),  satisfies  the 
consistency  condition  provided  that  the  original  difference  scheme,  Eq.  (80I, 
is  consistent.  Under  more  restrictive  assumptions,  they  were  able  to  show 
that  the  stability  of  the  ADI  scheme  follows  from  that  of  the  original 
difference  scheme,  and  also  that  the  final  solution  j?  ^approximates  ?r 
with  ay  error  no  worse  than  0 (At)" , and  consequently,  $ can  be  accepted 
as  ^ . Although  these  results  are  not  of  sufficient  generality  to  encompas.- 

the  Navier-Stokes  equations,  it  is  often  suggested  (e.g. , Ref.  67,  p.  .’15) 


,n+l 


67.  Richtmyer,  R.  D.  and  K.  W.  Morton:  Difference  Methods  for  Initial  Value 

Problems.  Second  Edition.  Interscience  ■>  New  York,  ft  i'J . 
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that  lit'  scheme  i:  st  th!  i . ■ j..-  ... . r ...  ttiv  us  more  m a : al  • 

those  for  w Lgorous  pi  . . ■ . 

the  procedure  as  outlined  here  is  appl : caul  • the  Mavier-Stoke.  aqu-oi  r 
is  adopted  as  a working  hypothesis,  support  x r th  .-  hyp  tt.ssis  is  > v 
by  the  favorable  result-  obtained  in  xctual  .•-om,-  utat  1 ons  ( ••••.•  Ref.  -3 

The  effort  involved  in  the  actual  programmi..,/  and  e - ;-a  a Rq 

di)  is  greatly  reduced . and  the  computer  storage  tvqui  • ..lem . . 

■ itractiog  Eq  <8l  from  Eq.  ( md  Eq.  (&  j m Sq.  Mat  E 

(8  ~ •)  ha\ • • the  following  sitn-li  fled  form: 

r(i^£)=B  ur-n 

s"(£'sr)  f) 

■ # ns  to  the  pi 

soluti  ms,  $ and  } , resi  cti  md  erve  to  stabilize  Eq.  (8l),  whici 

con!  tency  w the  go\  ai ng  equations.  r thi 
X .nenko  (Rex  65,  has  termed  the  Douglas -Gunn  procedure  the  method  o. 
stabili sing  correction  - . 

On  examination,  it  can  be  -eon  that  the  difference'  equation.  ( 
are  linear  in  the  *-level  quantities-.  At  tht  kth  step  in  the  prcwediu a . x 
are  m equations  at  each  grid.__pr.dnt  (xj,  x );  because  ox  the  s;  • • dif- 
ference operators  (6^  and  6^  ) tries e equation?  c attain  ti»e  depend et: • • a.  iasl*-. 
at  xk  and  at  each  of  the  two  adjacen*  grid  p ints  in  'he  x.  -dir- cti.  ..  fc  . 
sequently,  the  dix  ference  equations  mu;  r-  s iced  as  an  it'.ulioit 
should  be  recognized  that  upon  appli  ■ i ..  : a ..ucoe  -.  nun.:  - ; . . 

xR.  each  Of  r. ..  (8l.  84,  85)  generates  •*  block-1  . . 

....  - ' . • ' ' . ‘ ' 

Ai'ter  aj propria t.  treatment  of  c .xr.du:  uditi  ....  , ea.  n 

etfieientLy  using  a standard  block  eliminati.  r method  x-.ici  ax  • 
tact  rizat.ion  method  (:  ee  e . g . . R*.  f . 6l , . The  a tf  3d  ...  • 1 1 . • .. 
study  (Ref.  15)  is  clo.  ly  related  to  th.  matrix  fact  ri-uMun  tneM.. 

J 

the  elements  of  the  tridiagonal  matrix  .re  ^m  x u ibma- ri  o*  s r-  •!•••• 

1 n • lumber  of  1 


The  solution  procedure  for  a sint 1 

e t in.. 

step  is 

as  follows : 

(1  During  the  first  step  of  the 

ATI  j-roc.  iixrc. 

successive  Xj -direction  rowt 

of  grid 

implicit  ...  , . , 

The  implicit  systems  can  be  arranged 
olock-tridia.-owT  fora  v.l  solv  .1  as  1;.  ij  i ■ 

• there  are  • governing  w.iaticns.  1 i . r ii  . 

. yj  terns  have  (m  i . id  square  matrices  . th  vl  ,:i  mcu.tu  . 

-3)  fhe  second  and  third  steps  are  similar  to  , i exes  ’ tru.r  Eq.;. 

(84)  and  (05 ) are  applied  along  succes.  ive  x . nd  >:  direction 

rows  of  grid  point-  to  obtain  implicit  equo  Ions  tor  /**  and  $***. 
resi  ecti  . -fly . 

The  entire  solution  procedure  requires  cily  tv.  .'  levels  of  s* ... ■ ....  --  • ea.-h 
or  the  dependent  variables,  as  well  • i k-tridia*  ni  : 

i ■ . ...  im  1 is  required 

the  second  and  third  steps  in  the  p;  ;eedure.  when  Eqs.  (b4)  and  (to1  are 
...  I S Ref.  1 for  furthei  ■ , l . 


Further  Computational  Details 

Arti  ■: ' i . l .1  Visco.-ilr  T --c. 

:he  optic-:,  il  artificial  viscosi*.\  terms  which  wvr  tan-  , t uili> 

■ nee  equations  may  be  useful  in  predict  1 calculate  rr  t - stabilise  the 
metb-c  when  boundary  conditions  are  tr  a ted  in -ccuratel y,  when  coarse  mesh 
ut-acii.g  is  used , or  in  the  presence  01  discontinuities.  The  treatment  f 
artificial  viscosity  terms  used  here  (Ref.  i'5)  is  convt-M ■ v hu1  is  considered 
- r visio  . ■ 1 . ‘'nctliei  combination  cf  spatial  difference  formulas  and/or  arti- 
ficial viscosity  terms  may  eventually  prove  snperi-  however,  this  area 
requi  r e . f .r  • : ic.  s tudy . 

it.  i::  known  that  some  fin  it  ■ -difference  ap  pr  ox  j <t.a  t.  i on  c tc  fir.  r -r 
- - ’•  ition  errors  whi  ei  tab  Hi 

.ii  fferet.ee  . .item'-..  (Refs.  68  and  69).  It  seems  -.dvarr  . . 1 1 1 3 v f 

c in;-  • in.  vise. -us  flows,  to  use  a basic  iiffet  :iov  vo ...  . ...p 

ficial  viscosity  and  to  add  artificial  vJ  cosi’y  term.-'  it:  expi  i it  form  f r 
nal  use  w ne<  M tability.  * ..... 

central  Jifferences  for  ...  itial  derivative,  t.-  ..j;  -bit  from  41  ; < 

and  this  approach  has  Wen  taker. . Tl  • • d . -ri  >.  »r  1 fid  d vi  ■ • ,i*y  r. 


1 Roache,  P.  J.;  . t i ( 

Physics  Vol.  10.  197’,  r.  D-i. 
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suggested  by  a consideration  of  the  nondimensi<i:-'.l  model  eau-- 'ion 


u AAi  -1 
1 ax,  ' Re  dx,? 


(86) 


which  represents  a balance  between  advection  and  diffusion.  Roache  (Refs.  68 
and  69)  has  shown  by  a comparison  with  an  exact  solution  of  Eq.  (86)  that 
solutions  obtained  using  central  differences  for  the  advection  term  are  well 
behaved  provided  the  mesh  Reynold's  number  ReAx^  I ! yL  ^ but 

that  qualitative  inaccuracies  (associated  with  boundary  conditions  may  occur 
when  ReAx^  >2.  This  suggests  the  use  of  an  artificial  viscosity  to  ensure 
that  the  local  effective  mesh  Reynolds  number  is  no  greater  than  two.  Thus, 
Eq.  (86)  is  replaced  by 


From  Eq.  (88),  it  is  apparent  that  the  artificial  viscosity  can  be  made 
to  vanish  by  refining  the  mesh.  Since  the  artificial  viscosity  is  propor- 
tional to  Ax^,  solutions  will  have  first-order  formal  accuracy  if  Re^x j > 2 
but  second-order  accuracy  if  ReAx.<  2.  Strictly  speaking,  the  overall  method 
is  second-order  accurate  since  Re/^  -»o  as  the  mesh  is  refined.  It  should 
be  remembered,  however,  that  such  asymptotic  truncation  error  estimates  are 
meaningful  only  for  sufficiently  small  mesh  size;  whereas,  in  practical  cal- 
culations of  complex  flows,  mesh  resolution  capabilities  have  often,  out  of 
necessity,  been  strained.  One  virtue  of  the  present  formulation  is  that  by 
isolating  the  artificial  viscosity  terms  for  comparison  with  other  terms  in 
the  equations,  a nonrigorous  but  plausible  a posteriori  indication  of  the 
first-order  truncation  error  in  a computed  solution  is  available.  It  is,  of 
course,  obvious  that  Eq.  (86),  upon  which  the  artificial  viscosity  is  based, 
represents  a gross  simplification  of  the  Navier-Stokes  equations,  and  it  i- 
primarily  for  this  reason  that  the  present  formulation  of  artificial  ircosi^y 
terms  is  considered  provisional. 

Solution  of  the  Difference  Equations 


J 


It  has  been  shown  that  difference  equations  for  each  r-w  of  ►rrid  poin' 
in  the  three-dimensional  case  can  be  writ. ten  as  a block-tridiagonal  system 
having  (m  x m)  square  matrices  as  the  block  elements.  For  the  chemically 
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reacting  flow  considered  in  this  study  the  number  of  governing  equations  is 
as  least  seven  (m  £ 7).  Briley  and  McDonald  (Ref.  15)  have  shown  that  it 
is  possible  to  take  advantage  of  the  special  nature  of  the  coupling  in  the 
Navier-Stokes  equations  in  order  to  reduce  the  size  of  the  block  elements  in 
the  tridiagonal  system  which  must  be  solved  in  the  solution  procedure.  It  is 
only  necessary  to  solve  a block-tridiagonal  system  having  (m-2)x(m-2)  block 
elements  as  well  as  two  simple  tridiagonal  systems.  This  can  be  seen  by 
careful  examination  of  Eqs.  (1),  (4),  (5),  (30)  and  (31).  During  the  kth  step 
of  the  ADI  procedure  (k  = 1,  2,  3),  only  derivatives  with  respect  to  x^  and  t 
appear  in  the  implicit  difference  equations.  Therefore  only  the  x -direction 
momentum  equation  is  coupled  to  the  remaining  equations  (1),  (5),  (30)  and 
(31 )>  since  derivatives  of  the  other  two  velocity  components  with  respect  to 
x^  do  not  appear  in  those  equations. 

The  computational  effort  saved  by  solving  one  5x5  block-tridiagonal  and 
two  simple  tridiagonal  systems  rather  than  a 7x7  block-tridiagonal  system  can 
be  estimated  as  shown  in  Ref.  15.  If  there  are  N grid  points  along  a row 
being  treated  implicitly  and  l coupled  equations  at  each  grid  point,  then  the 
resulting  block-tridiagonal  system  requires  (3N-2)  (!?  + £)  arithmetic  opera- 
tions while  (5N-4)  operations  are  required  for  the  simple  tridiagonal  system 
(Ref.  6l) . Thus,  for  seven  coupled  equations  the  computational  effort  is 
approximately  llSON  and  for  five  coupled  equations  and  two  uncoupled  equations 
it  is  about  460N.  Therefore,  the  reduction  to  a system  of  five  coupled  equa- 
tions and  two  uncoupled  equations  during  each  ADI  step  affords  a considerable 
savings  in  computational  effort.  The  results  of  Ref.  15  indicate  that  the 
computational  effort  per  time  step  for  the  MINT  method  is  only  about  twice 
that  of  most  explicit  methods,  so  that  the  gain  in  efficiency  achieved  by 
taking  large  time  steps  is  real. 

Nonuniform  Grid  Transformation 


The  accuracy  of  solutions  computed  with  a given  number  of  grid  points 
can  often  be  improved  by  using  a nonuniform  grid  spacing  to  ensure  that  grid 
points  are  closely  spaced  in  regions  where  the  solution  varies  rapidly  and 
widely  spaced  elsewhere.  In  a typical  combustor  flow  situation,  large  gra- 
dients will  occur  near  the  combustor  injection  ports,  in  boundary  layers  on 
the  walls,  and  in  free  shear  layers  resulting  from  the  flow  separation  near 
injection  ports.  Transformations  for  the  latter  case  have  not  been  investi- 
gated in  the  present  study.  Boundary  layers  in  a turbulent  flow  are  treated 
using  the  wall  function  approximation,  although  occasionally  some  grid  trans- 
formation may  still  be  required  near  a wall. 


An  analytical  coordinate  transformation  has  been  devised  by  Roberts 
(Ref.  70)  which  is  a very  effective  means  of  introducing  a nonuniform  grid 
when  the  steep  gradients  occur  near  the  computational  boundarie:  . Suppose 
that  N grid  points  are  to  be  used  in  the  range  r X < X,,  and  that  steep 
gradients  are  anticipated  in  a region  of  thickness  6 (X.-X^)  near  X( . Then 
Roberts'  transformation  X^(X)  is  given  by 

XT(X)=N+(N-I)  m(^^)/ln(^)  (89) 

2 2 

where  a = X -X^,  b = a /(l-B),  and  c = X.  The  use  of  equally-spaced  points 
in  the  transformed  coordinate,  X..,,  ensure.:  an  adequate  resolution  of  both  the 
overall  region  X < X < X > and  the  subregion  X,  < X B(X. -X^ ) . Derivatives 
with  respect  to  the  physical  coordinate,  X,  are  obtained  from  the  following 
formulas : 


d . dxT  d (90) 

<?x  * dx  dxT 


The  use  of  three-point  difference  operators  for  X^  derivatives  in  Eqs.  (90-9.1 
produces  similar  operators  for  X derivatives.  These  X-derivativ-  operators 
can  be  computed  at  the  start  of  a calculation  and  stored,  along  with  the  • 
locations  of  grid  points.  No  further  consideration  of  the  transformation  is 
then  required. 


70.  Roberts.  G.  0.:  Computational  Meshes  for  Boundary  Layer  Problems. 

Proceedings  of  the  Second  International  Conference  on  Numerical  Methods 
in  Fluid  Dyanmics,  Springer-Verlag,  New  York.  1971.  p.  171. 


SECTION  IV 


KF  CULTS  AND  CONCLUSIONS 


T(.  evaluate  the  tnree-dimensional  computational  procedure  described 
in  the  : revious  sections,  comparisons  were  made  between  numerical  predic- 
tions obtained  fr  m the  MINT  code  and  unpublished  experimental  iata  tar.en 
in  a re  ■tan.n.ilar  research  combustor  by  the  iratt  Sc  Whitney  Aircraft  (i&WA 
C mb us t ion  Croup.  Toe  i&WA  results  consisted  of  temperature  measurements 
taken  with  a shielded  thermocouple  probe  and  emission  measurements  (unburne  : 
in  ir  rarbons,  nitric  oxide,  carbon  monoxide,  carbon  dioxide)  taken  using 
as  sampling  probes.  No  other  s\u table  measurements  were  available  for 
•mparison.  The  F&WA  research  c mbustor  is  a rectangular  duct  with  a cross 
section  1.5  by  3.0  in.  and  an  overall  length  of  about  10. i in.  The  flame 
h ’lder  employed  for  the  measurements  described  herein  was  a steel  plate 
( .25  in.  thick)  containing  an  array  of  eight  (3)  holes  as  shown  in  Pig.  3- 
The  holes  are  ' . 397  in.  in  diameter  with  0.75  in.  between  hole  centers. 

The  " &WA  combustion  experiment  employed  prevaporized  Jet-A  aircraft  fuel 
which  was  premixed  with  a preheated  air  stream.  The  nominal  fuel  composi- 
tion was  85.9  weight  percent  carbon  and  lL.l  weight  percent  hydrogen  wi‘ 
a molecular  weight  of  160.0;  therefore,  the  nominal  chemical  formula  for 
the  fuel  is  Cpp,l(L%2  38 • The  fuel  heat  of  combustion  is  approximately 
•t . 3x1  7 joule /kg  (18500  Btu/lbm),  and  a constant  fuel  specific  heat  (c, 
of  2510.0  .!oule/kg-°K  (0.6  Btu/lbm-°R)  was  used  in  the  present  calculations. 
The  air  preheat  temperature  measured  upstream  of  the  flame  holder  was 
about  756° K (900°F);  the  total  mass  flow  rate  through  the  combustor  was 
2.86xlC"2  kg/sec;  and  the  fuel  mass  fraction  was  0.0322. 

The  computational  region  boundaries  considered  in  the  present  calcula- 
tions are  indicated  by  dashed  lines  in  Fig.  3-  In  order  to  reduce  the 
number  of  grid  points  required  in  the  calculation  the  hole  pattern  was 
assumed  to  be  periodic  in  the  longitudinal  direction,  so  that  symmetry 
boundary  conditions  could  be  applied  and  only  one-half  of  the  inlet 
port  need  be  considered,  as  shown  in  Fig.  3-  The  temperature  measurements 
n the  combustor  indicate  that  this  is  a reasonable  approximation. 

The  computational  region  and  coordinate  system  are  illustrated  in 
Fig.  A.  Figure  5 shows  the  sections  (A,  B,  C)  for  which  profile  plots  of 
axial  velocity,  temperature,  and  nitric  oxide  are  presented.  All  predictions 
presented  were  ma  :e  using  a 17  by  9 by  17  grid  (26c  1 grid  points)  f r t.he 
■.  I , X2,  X3  directions,  respectively . A maximum  . i al  velocity  (Up)  of  lib. 7 
m see  (382.  o ft  sec ) was  required  to  match  the  e -perimental  mass  flow 
rate . The  referen’e  length  (h)  for  all  coordinates  is  C.Oloo'i  meters 
( .75  in.).  The  -alculatien  was  initiated  downstream  >f  the  flame  holder 
assuming  a mixture  inlet  temperature  of  75n°K  nn'  a pressure  of  1.002xJ.; 
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ia  (atmospheric  pressure).  Axial  velocity  profiles  are  riven  in  lips.  6, 
7,  and  8 at  sections  A,  B,  and  C,  respectively,  for  a series  of  axi' 1 
stations.  The  qualitative  characteristics  f the  et  expand i up  int  the 
combustor  seem  to  be  represented  quite  reasonably.  As  one  proceeds 
away  from  the  inlet,  the  jet  spreading  is  evident  and  reverse  fl  w is 
clearly  present  behind  the  inlet-plane  'wall  regions.  Due  t the  large 
inlet  velocity  and  the  relatively  small  number  of  grid  \ ints  em;  . yed 
in  the  calculation,  cell  Reynolds  numbers  in  the  axial  iirecti  n became 
large  so  that  significant  artificial  viscosity  was  required  in  the 
difference  equations  (see  SECTION  III).  The  presence  of  artificial 
viscosity  in  the  present  case  resuited  in  a 25  percent  increase  in  mass 
flux  between  the  inlet  and  computational  exit  plane  (x,  h.  , nondimen- 
sional).  This  problem  will  be  discussed  further  after  presental  n : 
the  remaining  results . 


The  equilibriur  hydr  -arbon  chemistry  assumption  was  not  expected 
to  be  valid  near  the  inlet  port  of  the  combustor  under  c ms 'deration. 
Therefore,  an  ad  hoc  ignition  delay  criterion  was  imposed  on  the  solution 
to  prevent  unrealistically  large  temperature  gradients  from  occurring 
near  the  inlet.  This  was  accomplished  in  the  present  numerical  procev . re 
by  specifying  the  pseudo-kinetic  rate  constant,  k;jC,  as  a function  of 
axial  distance  from  the  inlet  plane.  The  effect  of  the  specified  igniti  i 
delay  on  the  port-centerline  axial  temperature  variations  is  shewn  in 
Fig.  95  along  with  experimental  measurements  at  four  axial  stati  ns. 

The  gas  temperature  of  75b°K  (d  1 °]'\  originally  quoted  upstream  of 
the  flame  colder  is  shown  on  Fig.  -).  Subsequent  thermocouple  measurement r. 
with  similar  flame  holders  have  sh  wn  that  this  temperature  measurement 
is  probably  in  err  r and  should  be  a. proximately  1 30°K  (lh'0°F  , which 
is  consistent  with  the  measurement  at  - .1(  7.  Since  the  equilibrium 

temperatures  achieved  d vnstream  are  predicted  correctly  (Fig.  0,  it  seems 
i.x  ly  that  the  actual  temperature  far  upstream  of  the  flame  holder  is  75>  ' K. 
The  nigher  temperature  of  1 30° K believed  to  be  present  u stream  f the  flame 
holder  may  be  attributed  to  heat  transfer  through  the  flame  •.  L ter  fr  m the 
combustion  zone  to  the  gas  mixture  upstream  of  the  flame  ^ 

to  combustion  upstream  of  th<  f ame  holder.  Unfortunately,  tMi  Li  ra  i 
was  acquired  subsequent  to  completl  n of  tiie  present  numerical  -a Iculati. >ns , 
and  the  reported  predictions  are  thus  based  on  an  incorrect  inlet  temperature  . 

Proper  treatment  if  the  flame  holder  boundary  conditions  w i require 

additional  informat i m regarding  the  heat  transfer  or  surface  temperature  f 
the  flame  holler.  If  tnere  are  significant  energy  losses  thr  ugh  the  flame 
‘solder,  these  would  not  be  correctly  rc;  resented  by  the  adiabatic  wall 
conditions  actually  employed  in  the  calculations.  However,  the  adiabatic 
wall  conditions  are  consistent  xvith  tiie  inlet  temperature  of  7‘ "°K  used  in 
the  calculations.  One  further  consideration  affecting  the  temperature 
pre.iictions  is  the  error  in  mass  flux  attributable  ta  artifi  -ial  vise  sit y. 


Errors  in  mass  flux  associated  with  density  would  be  reflected  in  the 
temperature  through  the  equation  of  state.  Furthermore,  since  the  actual 
hydrocarbon  combustion  is  a none qui librium  process,  only  qualitative 
comparisons  are  warranted  between  the  predictions  and  measurements  where 
chemical  equilibrium  does  not  exit.  In  fact,  the  difference  between  compute : 
and  measured  axial  temperature  variation  in  Fig.  9 is  believed  to  be  due 
mainly  to  the  nonequilibrium  nature  of  the  combustion  chemistry,  and  thus 
accents  the  shortcomings  of  an  equilibrium  chemistry  model. 

Nondimens ional  temperature  profiles  are  presented  in  Figs.  1.  to  12 
at  sections  A,  B,  and  C,  respectively,  for  a series  of  axial  stations. 

The  presence  of  artificial  viscosity  (in  the  axial  direction)  in  both 
the  species  and  energy  equations  may  have  distorted  the  transverse 
temperature  profiles,  since  outside  the  Jet  boundaries  artificial  diffusion 
is  significant  compared  to  convection.  Representative  experimental 
temperature  measurements  and  the  corresponding  predictions  are  shown  in 
Figs.  13  to  15.  An  axial  translation  of  two  predicted  temperature  profiles 
(Figs.  l4  and  15 ) indicates  reasonable  qualitative  agreement  with  the 
measured  profiles.  This  translation  is  an  a posteriori  attempt  to 
compensate  for  the  incorrect  pseudo-kinetic  rate  constant  employed  in 
the  calculations.  The  qualitative  agreement  displayed  in  Figs.  14  and 
15  provides  some  confidence  in  the  numerical  procedure  and  justification 
for  further  improvement  of  the  physical  models.  Also,  the  distortion  iue 
t excessi  e artificial  viscosity  is  a problem  requiring  further  study 
before  quantitatively  meaningful  comparisons  can  be  made.  Both  Figs.  9 
and  r show  good  agreement  between  prediction  and  experiment  cnee 
chemical  equilibrium  has  been  achieved. 

Only  a limited  number  of  nitric  oxide  (NO)  concentration  measurements 
were  made  in  the  combustor  under  consideration.  Since  tiie  NO  concentration 
is  sensitive  to  the  temperature  history  of  a fluid  particle,  the  only 
comparison  presented  is  that  where  equilibrium  temperatures  have  been 
achieved.  The  predicted  NO  concentrations  (Fig.  16)  are  significantly 
lower  than  measured  values,  which  may  be  partly  attributed  to  the  simplify- 
ing assumptions  made  in  the  present  NO  chemistry  analysis  (SECTION  II)  and 
to  nonequilibrium  chemistry  effects  on  combustor  temperatures.  Furthermore, 
NO  concentrations  of  only  several  ppm  are  difficult  to  measure  accurately. 

The  present  numerical  results  demonstrate  both  the  basic  integrity 
of  the  computational  procedure  developed  under  the  present  study  an  i the 
capability  of  the  MINT  code  to  perform  calculati  ns  of  three-dimensional 
combusting  flows  with  recir  -illation.  The  need  for  improvement  of  the 
differencing  scheme  is  apparent  because  of  distortion  in  the  ; relictions 


due  to  excessive  artificial  viscosity.  Experience  at  UTRC  with  the  FREI 
computer  code  (Ref.  49)  has  shown  that  a two-equation  transport  model  of 
turbulence  offers  significant  advantages  over  a mixing  length,  eddy 
viscosity  turbulence  model  in  view  of  the  considerable  difficulty  of 
specifying  a mixing  length  distribution  in  a three-dimensional  recirculat- 
ing flow.  The  need  for  a simple  nonequilibrium  hydrocarbon  chemistry  model 
such  as  one  based  on  the  single  step  reaction  mechanism,  Eq.  (29) , is  also 
evident.  Work  on  alternate  difference  formulas  for  problems  with  large 
cell  Reynolds  numbers  and  on  a two-equation  turbulence  model  will  be 
carried  out  under  Phase  II  of  the  present  contract. 
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APPEND  O : JEKERALIZ!  Dili 


1 ' 1 Lneai  ■ : 

metric  coefficients  h^ , hp,  t«e  differential  elemeir  of*arc  length  da 
is  given  by 

(do  )?  - (h,  dx,)2  + (h2dx2)‘’  + (h3  dx3)2  J2  ) 

.e  vector  operat  i-oi.s  r -.•esrar,  ; r derivi  g • ne  g.  'verni  ,.r  different  Lai 
equations  in  the  present  coordinate  system  fron:  ue  vector  equations  1--  ) 
: elow.  Let  i^,  ig,  < . vectoi  in  the  co  >rd inat<  ] 

rections  x-^,  xp.  x respective!.  . ne:.  1 ::e  gradient  operator  is 

v _ _h d_  + j_e_  d_  ?3_  _d 

h|  dx i h2  dx2  + hj  dx3 


and  a general  vector  V is 


v-  'i  V'2V2  mivj 


The  divergence  of  a vec  n i_  give: 


h,  h2  h. 


3 [j77  <h2  h3  v|)  + ^.h3v2)+  (h,  h2v3)]  (<£) 


die  divergence  of  a synune  ric  -nsor  t with  components  t.^.  is  a vector  with 
the  following  component  s • of.  1J): 


(V  t ), 


[jx,  (h2hx'„'+  dx;(hihi'.2) 


(hlVl3)] 


1 h,h2h3  Ldx,  i"2"3V  + dx2^inj'l2]  ,tg 

dt,  + t , 3 dh,  _ t?2  dh2  1 33  dh3 

h|h3  dx3  h,h2  dx  i h , h 3 dx. 
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d 

, . 1 

d 

2 h,h2h3  1 

dx,  <h2h3fl2 

)+  u 

+ • > , ' ' ‘ . 

1 23  dti2 

, »IZ  dh2 

* " dhj 

in  dh. 

+ h2h3  d*  3 

h,h  dx. 

. i . 

(v  i),  = h-  ' . 1~7~ 

^ n i h2  h3  L(3x| 


1 13  ah3 
+ h|h3  dx i + 


1 23 
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(^2  ^ 3 ^ 13 ) 

+ dx2 (hi ^ 3 '23) 

s dh3 

tn 

3h, 

3 ^X2 

hl  h3 

ax  3 

d 

* 22 
h2  h3 


dht 
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The  foregoing  relationships  are  subst  ituted  into  the  eq  ia<  ions  ' L ) to  (5) 
to  obtain  the  governing  equations  in  the  generalized  coord i - ave  r,  c'em. 

The  axisymmetric  coordinates  are  obtained  by  taking  hg  = r,  where  r is  the 
radial  distance  from  the  axis  of  symmetry  and  xg  is  the  angular  coordinate 
(Fig.  1). 

An  analysis  for  determination  of  the  metric  coefficients  hj  and  h, 
based  on  the  two-dimensional  incompressible  potential  flow  in  a duct  has 
been  developed  by  Anderson  (Ref.  71'',  together  with  a computer  program  ‘ - 
implement  the  analysis.  The  method  utilizes  orthogonal  streamlines  and 
velocity  potential  lines  as  the  coordinate  lines  for  constant  x,  and  x,. 
respectively.  Therefore,  it  follows  from  potential  flow  theory  that  x,  i 

O p 

Xg  are  solutions  of  Laplace  s equation,  i.e. , y x,  y x,  = 0.  By  an 
appropriate  normalization  of  and  x^,  Anderson  i.  Ref . has  si  ■■ 

the  metric  coefficients  hj_  and  hg  can  be  made  to  satisfy  the  relati  .-.ship 
hf  = hg  = l/Uj , where  Uj  is  the  nondimens ion al  inviscid  velocity.  Anderson’s 
solution  procedure  is  based  on  a representation  of  the  boundaries  of  the 
inviscid  flow  region  by  a large  number  of  straight  line  segmen4  .• . thus 
forming  a many-sided  polygon,  .ho  Schwartz-Christ  ffei  nf  raal  trai  f r- 
mation  is  then  employed  to  transform  this  polygon  into  a semi-infinite 
plane,  where  the  potential  flow  solutii  ■ is  kn  wn.  ■ ral  Lve  n : ■ ri  :a3 

solution  for  the  transformed  location:,  of  4 he  various  nodal  poin4  (o  rners 
of  the  polygon),  followed  by  the  known  Lnvers<  rhwarl  - hrisl  ff<  L trai 
formation,  provides  the  inviscid  potential  flow  and  coordinate  rys4en  for 
the  desired  geometry.  Uniform-velocity  boundary  conditions  are  specified 
at  bhi  inlet  and  exi  pi  f th<  inviscid  flow  region.  f nec<  ar 

lesiri  . th<  1 . i . rid  fl  >w  regi  • • rxl  led  ips t ri  an  r dov  treai 

f th<  r<  loi  to  be  used  Ln  ; lie  visc<  us  sol  ition.  1 i s :omputei  pj  ran 
employed  bj  riley  and  Me  nald  • rZ)  to  obtai  a rth 
curvilinear  coordinate  system  f r e .-'rapid  at  ion  f 4 ::e-  t hree-dimenr  ioual  ;•  ..  - 
ic  flow  ii  rurved  passages.  ■ analysis  f ersoi  • f . q 
Imp!  ented  withi  - - frsu  rk  f this  si  < ■ ■ r : is  required. 


71 . Anders  . . ser 's  tai  al  for  a S ite-Dif f erence  aticn 

of  irbulen*  .'wirling,  'empress  ible  flow  i-  Ax  is  ,rno-  r i • .•  wil- 

r it  s and  .'L  ' >iled  '.-.alls.  AA"  - -7^-';0,  Vol.  . 1 - U. 

rile.,  ...  ' . and  • ".aid;  'umpi.4  a*  lor  of  -roe-  irm-.sior.nl  r- 

:len*  .'ul  sot.i<4  -Low  Ln  ' urved  ."assays.  sited  o>-  ese&rr 

4 er  epor‘  . - 'll1  * - . "are5  1 - 
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NOTE:  COMPUTATIONAL  BOUNDARIES  DENOTED  BY  DASHED  LINES 


Figure  3.  Experimental  flameholder  configuration  for  three  dimensional 
rectangular  combustor 
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Figure  8 Nondimensional  axial  velocity  profiles 


O PREDICTION 
A EXPERIMENT 


temperature  along  port  centerline 


R 1 1 32  6 


Figure  10.  Nondimensiona!  temperature  profiles 


Figure  12  Nondimensional  temperature  profiles 


SECTION  B X,  * 0.5 
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L = 0 01905  m 

O PREDICTION.  X3  = 0 219 
A EXPERIMENT,  X3  - 0 167 


TEMPERATURE  (NOIMDIMENSIONAL) 


Figure  13.  Comparison  between  predicted  and  experimental  temperature  profiles 
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Figure  16.  Comparison  between  predicted  and  experimental 
nitric  oxide  (NO)  concentration  profiles 
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